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Abstract 

Over  the  past  decade,  Gravity  Gradient  Instruments  (GGIs)  -  devices  which 
measure  the  spatial  derivatives  of  gravity,  have  improved  remarkably  in  accuracy  due 
to  the  development  and  rehnement  of  a  variety  of  accelerometer  technologies.  Some 
specialized  GGIs  are  currently  flown  on  aircraft  for  geological  purposes  in  the  mining 
industries  and,  as  such,  gravity  gradient  data  is  recorded  in  flight  and  detailed  gradient 
maps  are  created  after  post  mission  processing.  These  maps,  if  stored  in  a  database 
onboard  an  aircraft  and  combined  with  a  GGI,  form  the  basis  for  a  covert  navigation 
system  using  a  process  known  as  the  map  matching  method.  This  system,  if  it  could  be 
successfully  implemented,  would  be  completely  passive  -  impervious  to  conventional 
jamming  methods  and  relying  only  on  local  gravity  gradient  measurements  from  an 
onboard  sensor. 

This  paper  entails  an  investigation  into  the  feasibility  of  using  a  modern  GGI 
on  an  airborne  platform  for  covert  navigation  and  terrain  avoidance  by  examining 
GGI  signal  levels  in  different  flight  scenarios  (low,  medium,  and  high  altitudes  and 
velocities).  Previous  studies  using  gravity  gradiometers  have  been  accomplished  with 
promising  results  (some  theoretical  gradiometers  have  been  predicted  to  produce  GPS- 
like  navigation  accuracy).  However,  while  major  improvements  have  been  made  to 
current  airborne  gravity  gradient  instruments,  they  still  produce  noise  at  least  an 
order  of  magnitude  too  high  for  useful  aircraft  navigation  purposes.  This  research 
focuses  on  the  implementation  of  an  new  airborne  GGI,  currently  in  flight  test,  which 
has  demonstrated  approximately  an  order  of  magnitude  better  sensitivity  than  current 
airborne  GGIs.  To  demonstrate  whether  or  not  this  technology  is  currently  feasible, 
a  model  of  the  GGI  sensor  was  developed  to  investigate  signal  levels  at  representative 
flight  conditions.  Using  the  sensor  model,  representative  aircraft  trajectories  were 
flown  (simulated)  over  modeled  gravity  gradient  maps  to  determine  the  utility  of 
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flying  current  GGIs  in  the  roles  of  terrain  avoidance  and  navigation.  The  results  of 
the  GGI  simulations  at  different  altitudes,  velocities,  gravity  gradient  map  resolutions 
and  gradiometer  sensitivities  are  presented  and  discussed.  It  was  shown  that  the  map¬ 
matching  navigation  system  based  on  this  new  instrument  has  the  potential  to  provide 
a  marked  improvement  over  a  non-aided  INS  in  some  cases  but  was  limited  by  the 
drop  in  gravity  gradient  strength  at  higher  altitudes,  particularly  in  areas  of  smooth 
terrain.  It  was  originally  hypothesized  that  the  GGI  could  also  be  used  for  terrain 
avoidance  due  to  the  rapid  signal  change  as  rising  terrain  is  approached.  However,  GGI 
gradient  production  rate  and  bandwidth  limitations,  along  with  the  inverse  nature  of 
the  terrain  avoidance  problem,  rendered  GGI  aided  terrain  avoidance  unfeasible  for 
the  time  being. 
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An  Investigation  into  the  Feasibility  oe  using  a  Modern 
Gravity  Gradient  Instrument  for  Passive  Airgraft 
Navigation  and  Terrain  Avoidange 

I.  Introduction 

Background 

The  ability  to  precisely  navigate  is  a  critical  enabler  on  the  modern  battlefield. 
It  is  essential  to  mission  accomplishment  for  aircraft,  land  vehicles,  naval  vessels,  and 
even  personnel.  Perhaps  most  importantly,  in  a  military  utility  sense,  is  it  allows 
fast  moving  (i.e.  airborne)  platforms  to  find  their  way  to  and  place  weapons  on  a 
target  with  minimum  collateral  damage.  As  they  navigate,  most  modern  military 
aircraft  and  munitions  rely  on  some  form  of  a  Global  Navigation  Satellite  System 
(GNSS)  for  position  updates.  This  system,  using  a  set  of  signals  from  independent 
satellites  to  triangulate  position,  is  proven  and  provides  the  needed  accuracy  for  most 
mission  objectives.  However,  these  satellite  signals  can  be  denied  by  physical  blockage 
(i.e.  inside  a  cave  or  deep  underwater),  jamming  or  by  destruction  of  the  satellites 
broadcasting  them. 

There  is  a  significant  amount  of  research  into  methods  to  precisely  navigate  in 
a  GPS  denied  environment.  Some  of  these  include,  but  are  not  limited  to,  the  use 
of  pseudolites,  terrain  referenced  navigation  (TRN)  such  as  Sandia  Inertial  Terrain 
Aided  Navigation  (SITAN),  Terrain  Gontour  Matching  (TERGOM),  Terrain  Pro¬ 
file  Matching  (TERPROM),  image  based  navigation,  and  inertial  navigation  systems 
(INS)  which  can  be  provided  with  position  updates  from  the  aforementioned  naviga¬ 
tion  methods  to  correct  drift  errors  [1,2].  Another  method  of  aircraft  navigation  that 
has  been  given  relatively  little  attention  over  the  last  25  years  is  by  use  of  a  device 
known  as  a  gravity  gradient  instrument  (GGI). 

A  gravity  gradiometer  is  a  device  that  measures  spatial  derivatives  of  the  earth’s 
gravity  “acceleration”  vector.  These  spatial  changes  in  earth’s  gravity  are  caused 
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by  the  fact  that  the  earth  is  elliptical  (rather  than  perfectly  spherical),  is  spinning, 
has  varying  terrain  features,  and  mass  densities  which  are  not  constant.  They  are 
very  small  and  require  a  great  deal  of  sensor  accuracy  to  properly  measure.  Over 
the  past  two  decades,  gravity  gradiometers  have  been  carried  in  aircraft  and  used 
with  rapidly  increasing  success  for  geological  surveys.  The  speed  at  which  these 
aircraft  can  fly,  combined  with  improved  sampling  rates  and  noise  reduction  features  of 
modern  airborne  GGIs,  allow  surveyors  to  map  the  gravitational  gradients  caused  by 
terrain  and  subterranean  anomalies  much  faster  than  their  ground  based  counterparts 
(as  well  as  reaching  areas  otherwise  inaccessible  by  land).  Gravitational  gradient 
maps  are  recorded,  processed  and  used  post  flight  to  increase  understanding  of  the 
earth’s  gravitational  held  and  for  kimberlite,  oil  and  other  valuable  natural  resource 
detection  [3,4].  If  received  GGI  signals  were  able  to  be  correlated  to  an  existing  map 
generated  by  a  survey  (or  by  theory),  a  basis  for  a  covert  navigation  system  could 
be  formed,  similar  to  TRN,  but  requiring  no  external  emissions,  no  susceptibility 
to  adverse  weather  conditions,  and  a  signal  that  is,  by  today’s  standards,  virtually 
impossible  to  jam  -  requiring  terrain  to  be  moved  to  “fool”  the  sensor. 

Before  continuing,  “precision  navigation”  must  be  dehned  as  it  pertains  to  the 
scope  of  this  research.  The  term  precision  navigation  is  sometimes  loosely  thrown 
around  when  describing  the  accuracy  of  navigation  systems.  Genturies  ago,  precision 
navigation  was  a  matter  of  arriving  at  the  correct  continent.  During  WWII,  the 
Norden  bombsight  made  “precision”  high  altitude  bombing  a  reality  by  placing  bombs 
within  a  30m  circle  from  an  altitude  of  6km  (under  ideal  circumstances)  [5].  With 
the  advent  of  GPS  and,  more  recently,  differential  GPS,  precision  navigation  has, 
once  again,  been  redehned  with  navigation  errors  of  less  than  Im.  While  this  may 
seem  impressive  today,  suppose  in  the  future  that  one  wants  to  navigate  a  micro  UAV 
through  a  building  or  maybe  even  through  an  air  conditioning  vent!  Glearly,  Im  of 
error  could  be  unacceptably  large  for  that  application.  For  this  research,  precision 
navigation  refers  to  the  level  of  accuracy  attainable  by  GPS  or  GPS  aided  systems. 


2 


Gravitational  Gradient 


According  to  Wellenhoff  and  Moritz  [6],  the  gravitational  potential,  V,  of  a  point 
in  a  gravitational  field  is  defined  as  the  work  done  per  unit  mass  by  the  pull  of  gravity 
to  bring  a  body  from  inhnity  to  that  point.  It  is  a  scalar,  zero  order  tensor  function. 
From  Newtonian  potential  theory,  the  gravitation  potential  at  a  point  in  a  cartesian 
coordinate  system  {x,  y,  z),  due  to  an  attracting  mass  distribution  having  the  density 
function  p{x' ,y’ ,z')  and  volume  n',  is  given  by  the  following  volume  integral: 


V  = 


(1) 


where: 

r  =  \/{x  —  x'Y  +  {y  —  y'Y  +  {z  —  z'Y  and  represents  the  distance  between  the 
element  point  Q{x',y',z')  and  the  computation  point  P(x,y,z)  . 

G  is  Newton’s  gravitational  constant  and  is  6.6742  ■  10~ 

dv'  =  dx'dy'dz'  and  is  the  volume  element. 


The  gravitational  force  vector,  F,  is  the  gradient  of  the  gravitational  potential  and  is 
given  by: 


F  =  yv 


d\^  d\^  dV_ 
dx'  dy'  dz 


(2) 


The  gravitational  gradient  tensor,  Idj,  is  the  second-order  tensor  of  the  gravitational 
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potential  and  is  given  by: 


/  d'^V 
dx^ 

d^V 

dxdy 

d^V  \ 

dxdz 

(  V 

^  XX 
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The  gradient  is  symmetric  and  it’s  trace  satisfies  Poisson’s  equation:  =  —47rGp. 

When  the  density  at  the  computation  point  is  zero  (i.e.  free  air),  this  equation 
becomes  Laplace’s  equation.  Thus,  by  Laplace’s  equation,  which  states  that  the  trace 
of  the  tensor  must  sum  to  zero,  and  symmetry,  this  9  component  tensor  has  only  5 
independent  components. 

It  should  be  noted  that,  according  to  Equation  1,  the  gravitational  potential 
decreases  linearly  as  r  is  increased.  Consequently,  the  gravitational  force  and  gravi¬ 
tational  gradients  attenuate  as  a  function  of  and  r^,  respectively. 

Note  that  “gravitational”  phenomena  have  only  been  addressed  thus  far.  Grav¬ 
ity  is  a  more  familiar  term  and,  as  it  pertains  to  objects  on  the  earth’s  surface,  stems 
from  the  combination  of  the  gravitational  force  vector  and  centrifugal  force  due  to 
the  earth’s  rotation  (also  a  vector).  These  combined  forces,  acting  on  a  unit  mass, 
constitute  the  gravity  vector,  g.  “Gravity”  is  the  magnitude  of  vector  g  and  carries 
units  of  acceleration.  The  typical  value  of  this  acceleration  over  the  earth’s  surface  is 
the  familiar  9.8m/ Like  gravitational  potential,  gravity  potential,  W,  exists  and  is 
also  a  scalar,  zero  order  tensor  function.  It  is  simply  a  combination  of  gravitational 
potential,  V,  and  centrifugal  potential  <h: 


W  =  V  +  4> 


(4) 
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with 


g  =  Viy 


dW  dW  dW 
dx  ’’  dy  '  dz 


(5) 


In  geophysical  applications,  a  rotating  ellipsoid  of  revolution  is  used  to  approx¬ 
imate  the  earth  and  is  assumed  to  be  an  equipotential  surface  of  a  normal  gravity 
held  with  potential  U.  The  difference  between  the  actual  gravity  potential,  W,  and 
the  normal  potential,  U,  is  called  the  disturbing  potential,  T : 

T  =  W  -U  (6) 


with  the  gravitational  disturbance  gradients  dehned  as: 


/  d^T 
dx^ 

d^T 

dxdy 
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dxdz 
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where  T^x  is  the  change  of  gravity  in  the  x  direction  while  moving  a  known  distance 
in  the  x  direction  and  Txy  is  the  change  of  gravity  in  the  x  direction  while  moving 
a  known  a  known  distance  in  the  y  direction.  The  remaining  gradients  are  dehned 

similarly. 

The  gravitational  disturbance  gradient  tensor,  like  the  gravitational  gradient 
tensor,  satishes  Laplace’s  equation  in  free  air  and  is  symmetric,  thus  giving  it  5  inde¬ 
pendent  components  which  carry  units  of  1/s^.  Since  the  magnitude  of  the  gradients 
is  very  small,  units  of  1/ns^  are  commonly  used.  These  units,  known  as  Eotvos(i?o), 
were  named  after  19th  century  Hungarian  physicist  Baron  Roland  von  Eotvos  and  are 
not  recognized  by  the  SI  system  but  are  commonly  used  in  the  geophysics  commu¬ 
nity  [7].  To  add  physical  meaning  to  the  unit,  1  Eotvos  is  equivalent  to  the  gradient 
of  a  gravitational  held  produced  by  10  grains  of  sand  at  a  distance  of  1cm  [8].  Since 
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the  normal  gravity  potential  and  its  gradients  are  known  for  a  specified  ellipsoid  (i.e 
WGS84),  only  the  calculation  (or  estimation,  as  it  will  turn  out)  of  Tij  is  required. 
To  understand  how  Tij  manifests,  consider  its  source:  Assuming  the  earth  is  approx¬ 
imated  as  an  ellipsoid  of  revolution  with  a  smooth  surface  and  an  assumed  constant 
density,  variations  in  the  surface  (terrain)  and  density  contrasts  within  the  terrain 
and  the  ellipsoid  will  cause  variations  in  addition  to  the  nominal  potential  and,  in 
turn,  create  the  gravity  gradient  disturbance.  Some  of  the  methods  to  predict  these 
disturbance  gradients  will  be  investigated  in  Chapters  2  and  3.  It  should  be  noted 
that  a  GGI  measures  the  actual  gravity  gradients,  but  the  nominal  gradients,  Uij, 
are  known,  slowly  changing  (spatially),  and  generally  treated  as  a  bias.  Thus,  com¬ 
putation  of  the  gravitational  disturbance  gradient  tensor  is  the  more  urgent  focus  of 
current  research. 

A  Gravitational  Disturbance  Gradient 

To  gain  insight  into  the  gravitational  gradient  caused  by  a  mass  anomaly,  an 
example  using  a  simple  rectangular  prism  of  constant  density,  shown  in  Figure  (1),  is 
presented.  The  prism  is  defined  by  the  vertices  at  (xl,  yl,  zl),  (a;2,  yl,  zl),  {x2,  y2,  zl), 
and  {x2,  y2,  z2)  with  the  coordinate  system  used  having  axes  parallel  to  the  prism  sides 
and  the  origin  at  point  P.  Beginning  with  Equation  (1),  the  closed  form  solutions  for 
the  hve  gravitational  disturbance  gradients,  observed  at  point  P,  caused  by  the  prism 
can  be  found  [9]. 
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Figure  1:  Rectangular  Prism. 
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where: 


Ap  is  the  density  contrast  between  the  element  and  compntation  point  medinm 
r  =  ^/(x-  XiY  +  {y-  VjY  +  {z-  ZkY 

Note:  The  coordinate  system  used  here  and  throughout  this  study  is  the  North,  East, 
Down  (NED)  system  where  positive  y  corresponds  to  North,  positive  x  corresponds 
to  East,  and  positive  corresponds  to  Down. 


125m 

Figure  2:  Hypothetical  Prism  Orientation  and  Dimensions. 


Figure  2  shows  the  orientation  and  dimensions  of  the  prism,  or  hypothetical 
brick,  used  to  demonstrate  the  disturbance  gradient.  The  brick,  having  a  constant 
density  of  1.5p/cm^,  is  centered  on  a  250m  x  250m  grid.  The  gradients,  shown  in 
Figure  3,  were  calculated  on  a  plane  50m  above  the  center  of  the  brick  (z=-50m). 
Due  to  symmetry,  only  part  of  the  disturbance  gradient  tensor  is  shown.  The 
gradients  highlight  the  x-axis  (or  east-west)  edges  of  the  brick  by  measuring  the  east- 
west  changes  in  east-west  gravity.  Similarly,  the  Tyy  gradients  show  the  y-axis  (or 
north-south)  edges  of  the  brick  by  measuring  the  north-south  changes  in  north-south 
gravity.  T^z  highlights  the  overall  shape  of  the  anomaly  as  it  is  a  combination  of 
Txx  and  Tyy  with  a  sign  change.  T^z  and  Ty^  gradient  data  outlines  the  north-south 


Figure  3:  Example  Gravitational  Disturbance  Gradients. 

and  the  east-west  mass  anomaly  axes,  respectively.  They  also  help  to  highlight  the 
north-south  and  east- west  edges.  While  T^z,  Ty^  and  T^y  are  less  intuitive,  they 
contain  unique  information.  If  these  gradient  maps  were  stored  in  a  database  and  the 
gradients  were  able  to  be  accurately  measured  real-time  as  the  grid  were  traversed, 
there  is  enough  information  for  unique  determination  of  position  on  the  grid  based  on 
these  measurements.  This  is  the  fundamental  concept  behind  navigation  via  a  gravity 
gradiometer  based  map  matching  system. 

In  order  to  better  understand  the  frequency  content  of  the  gravitational  gradi¬ 
ents,  Figure  4  shows  the  spatial  frequency  spectrum  produced  by  the  brick’s  gradients. 
These  gradients  can  be  broken  into  spatial  frequencies  because  they  are  periodic  across 
position  in  space.  A  basic  understanding  of  an  anomaly’s  signal  structure  will  be  ben- 
ehcial  should  a  hlter  be  applied  to  a  real-world  gradiometer  signal.  The  plot  clearly 
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spatial  (cyc/m) 

Figure  4:  Example  Gravitational  Disturbance  Gradient  Spectrum. 

illustrates  that  most  of  the  signal  from  the  brick  is  in  the  0.001  to  0.01c|/c/m  spatial 
frequency  range  along  both  axes.  This  corresponds  to  wavelengths  of  approximately 
100  —  1000m.  Note  that  cyc/m  denotes  cycles  (or  periods)  per  meter  and  is  standard 
nomenclature  for  spatial  frequency. 


The  Gravity  Gradiometer 

Gonsider  a  proof  mass  attached  to  a  linear  spring  and  anchored  inside  a  housing 
in  a  reference  frame  which  is  free  from  a  gravitational  held  (see  Figure  5).  In  this  case, 
Newton’s  Second  Law  is  simply:  mx  =  F,  where  m  is  the  mass  of  the  proof  mass, 
X  is  acceleration  along  the  x  axis  and  F  is  the  force  applied  to  the  housing.  When 
a  specihc  force  acts  on  the  housing,  it  will  accelerate  with  constant  acceleration,  a, 
with  respect  to  the  given  reference  frame.  This  will  cause  the  proof  mass  to  move 
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Figure  5:  Accelerometer. 


relative  to  the  housing,  the  spring  to  compress,  and  the  resultant  spring  force,  fg, 

to  act  on  the  mass.  Let  Xpm  denote  the  position  of  the  proof  mass  relative  to  the 

housing  and  Xh  denote  the  position  of  the  housing  in  the  external  reference  frame. 

Now,  the  position  of  the  proof  mass  in  the  external  reference  frame  is  x  =  +  Xprn 

and  it’s  acceleration  is  x  =  x/^  +  Xpm-  The  spring  force,  given  by  Hooke’s  law,  is: 

fs  =  —kxpm  =  mxpm,  where  k  is  the  spring  constant.  Therefore,  by  Newton’s  Second 

k 

Law  of  Motion,  the  equation  of  motion  for  the  proof  mass  is:  Xp^  H - =  —a. 

m 

Assuming  initial  conditions  of  Xpmit  =  0)  =  0  and  Xpmit  =  0)  =  0,  the  solution  to  the 
proof  mass’s  differential  equation  is: 


(13) 


Thus,  the  position  of  the  proof  mass  relative  to  the  housing  is  proportional  to  the 
applied  acceleration  (with  proportionality  constant  ^).  If  or  x  can  be  measured,  a 
can  be  found  in  which  case  this  device  is  now  an  accelerometer.  While  this  constitutes 
a  simple  example  of  an  accelerometer,  devices  in  use  today  are  based  on  the  same 
fundamental  principles  (i.e.  somehow  measuring  the  relative  motion  of  a  proof  mass 
to  solve  for  acceleration).  Now,  suppose  this  simple  accelerometer  is  placed  in  an 
area  where  it  is  under  the  influence  of  a  gravitational  field,  but  no  specific  forces  act 
on  the  housing  (e.g.  freefall).  Newton’s  Second  Law  becomes  (assuming  rUg  =  mj. 
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that  is,  inertial  and  gravitational  mass  are  the  same):  x  =  a  +  g,  where  g  is  the 
gravitational  acceleration  vector.  This  gravitation  vector  will  act  to  accelerate  not 
only  the  housing  itself,  but  also  the  proof  mass  and  the  spring.  By  Newton’s  Law  of 
Gravitation,  the  gravitational  acceleration  of  the  housing  and  all  of  its  components  is 
the  same  (assuming  the  gravitational  acceleration,  g^  is  constant  over  the  housing). 
Now  the  motion  of  the  proof  mass  in  the  external  frame  becomes:  x  =  g.  Likewise, 
Xh  =  g.  Thus, 

X  =  Xh  +  Xp^  ^  g  =  g  +  Xpm  ^  =  0  ^  Xpm{t)  =  0  (14) 

and  there  is  no  motion  of  the  proof  mass  relative  to  the  box.  In  other  words,  the 
accelerometer  is  accelerating  in  a  gravitational  held,  but  measures  no  acceleration! 
That  is,  an  accelerometer  does  not  directly  sense  the  presence  of  a  gravitational  held, 
only  specihc  forces  resulting  from  applied,  action  or  contact  forces.  To  reiterate, 
accelerometers  do  not  sense  gravitational  acceleration.  They  will,  however,  sense 
reactions  from  gravitational  held  forces.  For  example,  if  an  accelerometer  oriented 
along  the  “down”  axis  in  a  NED  reference  frame  were  placed  on  the  surface  of  the 
earth,  it  would  sense  the  reaction  to  the  gravitational  force  provided  by  the  earth’s 
surface.  Finally,  in  accordance  with  the  principle  of  equivalence,  the  accelerometer 
cannot  distinguish  whether  this  reaction  is  a  result  of  gravitation,  rotation,  or  an 
applied  force.  This  trait  is  the  key  behind  the  concept  of  the  gravity  gradiometer  [10]. 

2  Accelerometers  at  a  known  distance  apart 

A2 

L - 1 

Figure  6:  A  Simple  Gravity  Gradiometer. 

Suppose  that  an  accelerometer  is  used  in  an  attempt  to  measure  gravity  reac¬ 
tion  forces.  This  device  represents  an  gravimeter  -  a  single  accelerometer  oriented 
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to  measure  gravity  along  some  axis  of  interest.  While  gravimeters  are  very  success¬ 
ful  in  measuring  gravity  information  on  fixed  platforms  (such  as  the  ground),  their 
accuracy  plummets  when  used  on  a  moving  base  platform  -  particularly  in  the  air¬ 
borne  environment.  Aircraft  vibrations  from  engines,  pumps,  etc.,  combined  with 
accelerations  from  turbulence,  engine  thrust  changes,  and  maneuvers  to  render  ac¬ 
curate  gravity  measurements  difficult  due  to  the  single  accelerometer’s  inability  to 
distinguish  between  inertial  and  gravitational  acceleration. 

Now  suppose  that  a  pair  of  accelerometers  are  mounted  in-line  along  some  axis 
at  some  known  distance  apart,  as  in  Figure  6.  If  the  accelerometer  readings  are 
differenced  and  then  divided  by  the  length  between  them,  a  gravity  gradient  has  been 
measured: 


gravity  gradient  = - - -  (15) 

Jj 

By  measuring  the  gravity  gradient,  host  vehicle  accelerations  of  the  first  order  are 
intrinsically  rejected,  thus  leaving  only  the  differential  acceleration  of  the  earth’s 
gravity  held  over  some  unit  distance  [11].  It  should  be  noted  that  the  distance  between 
the  accelerometers  is  critical  for  gradiometer  performance.  If  the  distance  is  too 
large,  the  host  vehicle  accelerations  sensed  by  each  accelerometer  may  be  substantially 
different  and  thus  difficult  to  difference  out.  If  the  distance  is  too  small,  gradiometer 
sensitivity  will  be  compromised.  In  reality,  accelerometer  misalignment,  scale  factor 
differences  and  other  noise  sources  can  corrupt  the  gradient  measurement.  In  order 
to  measure  the  full  tensor  of  gravity  gradients,  a  minimum  of  three  accelerometer 
pairs  are  oriented  along  three  axes.  It  should  also  be  noted  that  the  terms  “gravity 
gradient  instrument” ,  “GGI” ,  and  “gravity  gradiometer”  are  used  interchangeably  in 
this  research. 
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An  Inverse  Problem 


Suppose  that  the  shape  of  the  earth  and  the  density  variations  within  it  are 
exactly  known.  With  this  information,  a  unique  value  of  the  actual  potential,  W,  can 
be  found.  That  is,  the  determination  of  the  actual  potential  is  a  “direct”  problem. 
Now  suppose  that  the  actual  potential  is  known  (perhaps  by  GGI  measurements) 
and  the  shapes  and  density  variations  that  caused  it  are  to  be  calculated  (i.e  the 
“inverse”  problem).  It  is,  in  fact,  impossible  to  uniquely  solve  for  these  potential 
generating  masses  without  additional  information.  There  are  an  inhnite  number  of 
possible  combinations  of  mass  location  and  density  variations  that  could  create  a 
certain  value  for  the  potential.  To  determine  the  solution  (or  to  better  estimate  it), 
additional  information  must  be  provided.  This  phenomenon  may  hamper  GGI-based 
terrain  avoidance  performance. 

Problem  Statement 

While  a  limited  number  of  navigation  performance  studies  using  information 
from  a  gravity  gradio meter  have  been  accomplished  in  the  past,  the  high  error  level 
associated  with  using  this  sensor  on  an  aircraft  rendered  successful  navigation  mainly 
a  function  of  the  assumptions  regarding  the  performance  of  future  gradiometers.  Be¬ 
cause  of  these  errors  and  the  inverse  problem,  relatively  little  research  into  GGI  based 
navigation  and  terrain  avoidance  is  available  in  open  literature.  Since  airborne  grav¬ 
ity  gradiometers  have  seen  remarkable  improvements  over  the  past  decade  [12],  this 
research  aims  to  investigate  the  feasibility  of  using  a  modern  GGI,  which  must  pro¬ 
vide  real-time  gradient  measurements,  in  the  role  of  passive  navigation  and  terrain 
avoidance  with  emphasis  on  military  type  flight  environments  by  rigorously  examining 
simulated  GGI  signals  at  a  variety  of  representative  flight  conditions  and  comparing 
them  to  those  proven  in  previous  works  to  yield  navigation  performance  improve¬ 
ments.  The  signals  will  be  analyzed  with  respect  to  their  ability  to  be  matched  to 
a  map.  It  should  be  noted  that  most  of  the  representative  flight  conditions  to  be 
examined  have  not  been  previously  investigated  in  open  literature. 
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This  research  is  part  of  the  hrst  phase  of  a  three-phase  plan.  Phase  one  is  a  fea¬ 
sibility  study  where  a  GGI  will  be  selected  and  its  model  developed  and  implemented 
to  investigate  signal  levels  at  representative  flight  conditions.  Phase  two  involves 
procuring  a  GGI  and  integrating  it  as  a  sensor  in  a  navigation  system.  This  phase 
will  include  sensor  model  rehnements  and  validation  and  will  culminate  with  a  navi¬ 
gation  demonstration  on  a  land-based  vehicle.  The  hnal  phase  involves  flight-testing 
of  the  navigation  system  in  an  aircraft  to  demonstrate  military  utility  and  validate 
the  overall  modeling  effort. 

Feasibility 

First  and  foremost,  the  dehnition  of  feasibility  within  the  scope  of  this  study 
must  be  dehned.  The  following  stipulations  will  apply: 

•  Navigation  will  be  performed  onboard  an  airborne  aircraft. 

•  The  gravity  gradiometer  will  be  in  size  and  weigh  450kg,  maximum. 

•  “Modern  GGI”  is  dehned  as  a  gravity  gradiometer  projected  to  be  available 
within  the  next  10  years. 

•  Passive  Navigation  will  be  based  on  a  map-matching  method  and  performance 
improvements,  if  any,  will  be  measured  against  unaided  Nav  Grade  IMUs. 

•  Terrain  avoidance  performance  will  be  based  on  the  GGPs  ability  generate  a 
signal  that  is  useable  to  predict  terrain  impact  is  imminent  within  1.5  seconds 

[13], 

Research  Objectives 

With  feasibility  dehned,  the  research  objectives  for  this  study  are  presented. 
It  should  be  noted  that  sensor  cost  was  not  considered  for  this  study.  Additionally, 
further  assumptions  and  limitations  will  be  addressed  in  subsequent  chapters. 
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•  Investigate  modern  gravity  gradiometers  available  for  introduction  into  an  air¬ 
craft. 

•  Select  gradiometers,  based  on  noise,  bandwidth,  sampling  rate,  size  and  weight, 
most  suitable  for  aircraft  navigation  and  terrain  avoidance. 

•  Develop  models  of  the  GGI  sensors  deemed  most  usable  for  aircraft  navigation 
and  terrain  avoidance. 

•  Generate  gravity  gradient  maps  that  represent  realistic  values  produced  by  the 
earth. 

•  Examine  simulated  GGI  signal  variations  in  response  to  factors  including  alti¬ 
tude,  airspeed,  terrain  variation,  and  formation  effects. 

•  Attempt  to  classify  signal  threshold  levels  for  useful  terrain  avoidance  and  nav¬ 
igation  via  map-matching. 

•  Determine  if  the  selected  GGI  meets  signal  threshold  requirements. 

•  Recommend  needed  gradiometer  improvements,  if  any,  and  appropriate  ways 
to  integrate  the  GGI  signal  into  navigation  (with  emphasis  on  map-matching 
methods)  and  terrain  avoidance  systems. 

Preview 

This  thesis  is  divided  into  four  subsequent  chapters.  Ghapter  II  presents  the 
literature  review  for  this  research.  Divided  into  three  parts,  it  encompasses  the  history 
of  gravity  gradiometry,  a  review  of  modern  airborne  gravity  gradiometer  technology, 
and  a  review  of  previous  gravity  gradiometer  based  navigation  and  terrain  avoidance 
research.  Ghapter  III  highlights  the  problem  setup  and  methodology.  It  covers  how 
gravity  gradient  maps  were  constructed,  how  the  GGI  was  modeled,  and  the  tests 
that  were  executed  in  order  to  determine  navigation  and  terrain  avoidance  feasibility. 
Ghapter  IV  provides  the  results  and  analysis  and  serves  to  report  the  Endings  from 
the  feasibility  study.  Ghapter  V  is  a  closing  discussion  that  will  conclude  the  thesis 
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with  significant  contributions  and  insights.  Also,  some  challenges  and  future  research 
recommendations  will  be  discussed. 
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II.  Literature  Review 


History  of  Gravity  Gradiometry 

Gravity  gradiometry  began  in  1890  when  Baron  Roland  von  Eotvos,  a  Hungar¬ 
ian  physicist,  developed  an  instrument  known  as  the  torsion  balance  to  measure  the 
minute  variations  in  gravity  over  a  short  distance.  The  torsion  balance  was  made  of 
a  metal  beam,  suspended  by  a  wire,  with  weights  at  each  end  (similar  to  a  dumb¬ 
bell).  If  gravity  varied  with  position  along  the  axis  the  weights  were  placed  on,  the 
force  exerted  on  each  weight  would  be  different,  thus  causing  a  rotational  force  on  the 
beam  and  in  turn  causing  the  wire  to  twist.  Eotvos  measured  the  amount  of  twist  to 


Figure  7:  Eotvos’  Torsion  Balance  -  The  First  Gravity  Gradio meter 

In  1901,  Hugo  de  Boeckh,  head  of  the  Hungarian  geologic  survey,  convinced 
Eotvos  to  test  the  real-world  usefulness  of  the  torsion  balance.  The  device  was  used 
to  map  the  shape  of  a  frozen  lake  basin,  which  was  already  well  known  from  previous 
summertime  measurements  made  from  a  line  and  sinker.  The  test  was  a  resounding 
success  -  the  contour  map  generated  via  the  torsion  balance  matched  the  previously 
made  maps.  Eotvos  and  Boeckh  then  completed  more  difficult  geological  surveys  in 
the  region.  Word  of  Eotvos’  torsion  balance  success  quickly  spread  to  oil  prospectors 
around  the  world  -  gravity  gradiometry  had  officially  been  born. 

After  World  War  I,  American  geologists  used  the  torsion  balance  in  a  attempt  to 
find  salt  domes  -  mushroom  shaped  underground  geologic  structures  that  often  have 
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oil  and  gas  deposits  along  their  sides.  Since  salt  is  less  dense  than  most  rock,  it  exerts 
a  weaker  gravitational  force  relative  to  the  earth  surrounding  it.  As  such,  gravity 
gradients  can  highlight  a  buried  salt  dome.  In  1924,  geologists  from  the  Amerada 
Hess  Corporation  struck  paydirt  by  hnding  a  hidden  salt  dome  via  measurements 
made  from  a  torsion  balance.  By  1935,  the  use  of  gravity  gradiometry  for  subsurface 
surveys  was  routine  -  particularly  in  the  oil  business  [14]. 

The  early  success  of  the  torsion  balance,  however,  did  not  secure  its  long  term 
use.  The  instrument  was  fairly  difficult  to  use  in  the  held.  In  order  to  make  a  reliable 
measurement,  geologists  had  to  hrst  clear  a  100  meter  long  swath  in  eight  directions 
(star  pattern)  from  the  location  of  the  torsion  balance  to  prevent  the  mass  of  trees  and 
rocks  from  corrupting  measurements.  Additionally,  a  small  building  had  to  be  erected 
in  order  to  protect  the  instrument  from  wind  and  temperature  changes.  To  get  an 
idea  of  the  sensitivity  of  the  torsion  balance,  measurements  could  be  corrupted  by  the 
large  belt  buckles  often  worn  by  geologists!  To  compound  the  problem,  gradiometer 
data  was  often  misinterpreted  which  led  to  false  survey  results.  These  issues  led  to  the 
boom  in  the  use  of  gravimeters,  devices  which  measure  gravity  rather  than  the  change 
of  gravity  per  unit  distance,  for  surveys.  Gravimeters  are  inherently  less  sensitive 
than  gradiometers  and  thus  did  not  require  extensive  measurement  site  preparation. 
Furthermore,  the  data  from  gravimeters  was  easier  to  interpret.  This  led  to  increased 
investments  in  gravimeters  and  by  the  1950s,  gravimeters  had  replaced  gradiometers 
in  most  gravity  held  measurement  applications  [4].  For  the  time  being,  the  gravity 
gradiometer  was  gone  but  certainly  not  forgotten. 

In  the  1970s,  both  US  and  Russian  navies  realized  that  the  accuracy  of  sub¬ 
marine  launched  ballistic  missiles  (SLBMs)  depended  greatly  upon  precise  knowledge 
of  gravity  at  the  time  of  missile  launch.  Since  gravimeter  measurements  plummet  in 
accuracy  on  moving  platforms,  a  new  wave  of  research  into  gravity  gradiometers  was 
launched.  [4]  Around  the  same  time,  the  US  Air  Force  had  abandoned  gravimeter  sys¬ 
tems  for  airborne  surveys  due  to  the  fact  that  kinematic  accelerations  overwhelmed 
the  anomalous  gravimetric  signal  on  the  aircraft  in  hight.  By  the  early  1980s,  Bell 
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Aerospace  Textron  had  successfully  developed  a  moving  base  full  tensor  gravity  gra- 
diometer.  The  instrument,  developed  by  Ernest  Metzger,  was  selected  by  the  Navy 
(with  over  400  million  USD  in  development  costs)  for  gravity  compensation  require¬ 
ments  of  its  trident  submarine  inertial  navigation  systems  and  by  the  Air  Force  Geo¬ 
physics  Laboratory  (AFGL)  for  it’s  region  airborne  gravity  survey  system.  In  the  mid 
1980s,  part  of  this  technology  was  declassihed  and  eventually  used  in  the  1987  Defense 
Mapping  Agency  (DMA)  funded  flight  test  of  the  Gravity  Gradiometer  Survey  Sys¬ 
tem  (GGSS).  This  test,  accomplished  by  AFGL,  constituted  the  first  airborne  gravity 
gradiometer  survey  published  in  open  literature.  The  GGSS  consisted  of  Bell/Textron 
(now  owned  by  Lockheed  Martin)  gradiometers  which  were  installed  in  the  back  of  a 
van  along  with  other  support  equipment.  Since  much  of  the  GGSS  was  hardwired  into 
the  van,  the  van  was  simply  loaded  into  a  G-130  Hercules  aircraft  for  flight  test.  The 
survey  was  flown  over  southwestern  Oklahoma  and  northern  Texas  and,  while  high 
in  noise  (~  AOEo/y/H^),  was  able  to  measure  low  frequency  effects  corresponding  to 
subterranean  anomalies  in  the  area  [15].  In  all,  the  GGSS  represented  an  outstanding 
achievement  that  sparked  a  fury  of  airborne  gradiometer  development.  The  speed  at 
which  aircraft  can  fly,  as  well  as  the  ability  to  access  remote  areas  of  land,  made  much 
larger  and  quicker  surveys  a  reality.  Many  oil  and  other  valuable  natural  resource 
mining  industries  had  renewed  interest  in  gravity  gradiometers. 

Today,  airborne  gravity  gradiometers  are  used  mainly  for  geological  surveys  in 
the  hunt  for  valuable  natural  resources.  Gompanies  such  as  Bell  Geospace,  ARKex, 
Gedex  and  Fugro  provide  airborne  gravity  gradient  surveys  to  customers  who  desire 
such  data.  Also,  geophysicists  use  them  to  better  understand  our  planet’s  gravity 
held  and  overall  structure.  For  this  role,  gradiometers  have  been  installed  and  used 
with  success  on  ships,  aircraft  and  satellites  [4,16-18]. 
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Airborne  Gravity  Gradiometer  Technology 


Figure  8:  Rotating  Accelerometer  Schematic,  taken  from  [3]. 

Rotating  Accelerometer  Gravity  Gradiometer.  The  rotating  accelerometer 
Gravity  Gradient  Instruments  (GGIs)  are  based  on  the  Bell  Aerospace/  Textron  de¬ 
sign  (now  owned  by  Lockheed-Martin)  and  use,  at  a  minimum,  2  pairs  of  conventional 
accelerometers  monnted  opposite  of  each  other  on  a  spinning  disk  to  measure  gradi¬ 
ents  in  the  plane  of  rotation  (i.e.  normal  to  the  axis  of  spin)  as  shown  in  Figure  8 
Each  accelerometer  consists  of  a  mass  which  is  pivoted  (i.e.  a  pendulum)  and  a  sensor 
that  measures  the  offset  position  of  the  pendnlnm  along  its  path  of  travel.  Included 
within  each  accelerometer  is  electronic  circnitry  that  restores  the  pendnlnm  to  its  base 
position  throngh  the  use  of  electromagnets  and  constrains  the  pendnlnm  to  minimal 
movements  along  the  inpnt  axis  of  the  accelerometer  [19].  This  applied  electric  signal 
represents  the  ontpnt  of  the  accelerometer  and  serves  as  a  measnre  of  the  acceleration 
of  the  pendnlnm  bronght  on  by  any  forces  applied  to  the  accelerometer.  The  mea- 
snrements  from  each  pair  of  accelerometers  can  be  resolved  into  two  gradients  in  the 
plane  of  the  rotating  disc  by  acconnting  for  the  distance  between  each  accelerometer, 
the  rate  at  which  the  disc  is  spinning,  and  the  difference  in  the  measnred  accelerations 
between  each  pair.  In  order  to  obtain  a  fnll  gravity  tensor  (5  independent  gradients), 
three  rotating  discs  must  be  used  since  each  disk  can  only  measnre  2  components  of 
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the  tensor.  The  tensor  components  measnred  in  the  external  coordinate  axis  are  then 
fonnd  by  nsing  the  appropriate  linear  combination  of  the  six  GGI  ontpnts. 

The  three  rotating  accelerometer  GGI  ontpnts  are  given  as  follows: 

xz  plane  :  {Ai  +  A2)  —  (vda  +  A4)  =  2RdsinVtt{Tzz  —  T^x)  +  ^RdTxzCos2Vtt,  (16) 

yz  plane  :  {Ai  +  A2)  —  (vds  +  A4)  =  2RdSinVtt{Tzz  —  Tyy)  +  4:RdTyzCos2Qt,  (17) 

xy  plane  :  {Ai  +  A2)  —  (A3  +  A4)  =  2RdSinVtt{Tyy  —  Txx)  +  4:RdTxyCos2Q,t  (18) 

More  specihcally,  the  derivation  of  the  measnred  ontpnt  from  a  single  disk  is 
presented:  From  the  geometry  of  the  disk  and  placement  of  the  accelerometers  shown 

SPIN 


Fignre  9:  Single  Rotating  Accelerometer  Disk 

in  Figure  9,  an  equation  for  the  acceleration  measured  by  accelerometer  Ai  can  be 
derived  [11]: 


Ai  =  {ay  +  TyxRdCosflt  +  TyyRdSinQt)cosflt— 

{ax  +  TxxRdCosflt  +  TxyRdSinQt)sinQt  (19) 

where  ay  is  the  gravitational  held  induced  acceleration  at  the  center  in  the  y  direction, 
ax  is  the  gravitational  held  induced  acceleration  at  the  center  in  the  x  direction,  Rd 
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is  the  distance  from  the  center  of  the  disk  to  the  accelerometer,  and  hi  is  the  angular 
velocity  of  the  disk  about  its  spin  axis. 

Expanding  Equation  19  gives: 


=  {dyCosVtt  +  TyxRdCos^Vtt  +  TyyRdsinVtt  cosilt)  — 

{axsinilt  +  TxxRdCosilt  sinilt  +  T^yRdSin^ilt)  (20) 


Recalling  the  following  trigonometric  identities:  cos^il  =  \  +  ^cos2r2,  siri^il 
^cos2fl,  and  siniltcosilt  =  \sin2ilt  with  Tyx  =  T^y,  Equation  20  gives: 


Al  =  ttyCOSilt 


1  1  Rd 

axsinflt  +  TxyRd{-  +  -cos2flt)  +  Tyy—sin2ilt 


2 

Txx^sin2Qt  -  TxyRdi]:  -  ]rCos2Qt) 


'2 
Rd 
~2 


1 

2 


(21) 


Combining  like  terms  of  Equation  21  gives: 

R 

Al  =  ttyCosilt  —  OxSinilt  +  TxyRdCos2Qt  +  sin2VLt{Tyy  —  Txx)  (22) 

Since  the  opposing  accelerometer  {A2  in  this  case)  is  always  tt  radians  away  from 
Al,  the  acceleration  measured  by  A2  can  be  derived  by  replacing  ilt  with  fit  +  tt  in 
Equations  (l)-(4): 

R 

A2  =  —dycosilt  +  axsinilt  +  TxyRdCos2flt  +  —sin2ilt{Tyy  —  Txx)  (23) 


Summing  Equation  22  and  Equation  23  gives: 


Al  +  A2  =  2TxyRdCos2ilt  +  Rdsin2flt(Tyy  —  Txx)  (24) 

Replacing  ilt  with  flt+|  and  ilt+^  in  Equations  19-23  gives  the  relationship  between 
R3  and  A^: 

A^  +  A4  =  —2TxyRdCos2ilt  —  RdSin2ilt{Tyy  —  Txx)  (25) 
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Subtracting  Equation  25  from  Equation  24  gives  the  basic  element  of  measurement 
from  a  rotating  disc  gravity  gradient  instrument: 

{Ai  +  A2)  —  {Aj,  +  A4)  =  2RdsinVtt{Tyy  —  T^x)  +  ^RdTxyCos2VLt  (26) 


This  combination  signal  (Equation  26),  called  a  bandpass  signal,  is  normally 
low-pass  hltered  and  digitized,  and  then  demodulated  at  sin2flt  and  cos2Qt  to  obtain 
Txy  and  (Tyy  —  Txx)-  Also  note  that  if  the  accelerometers  are  perfectly  aligned,  scale 
factor  balanced,  and  linear,  no  angular  or  wheel  acceleration  terms  appear  in  Equa¬ 
tion  26.  Additionally,  any  residual  linear  acceleration  sensitivity  will  be  modulated 
at  n  and  will  not  appear  after  the  demodulation  at  2VL.  In  essence,  the  perfect  rotat¬ 
ing  accelerometer  gradiometer,  if  mounted  on  a  stabilized  platform,  is  not  sensitive 
to  vehicle  accelerations  to  the  hrst  order  [11].  However,  sensor  misalignment,  scale 
factor  differences  of  each  accelerometer,  and  other  real-world  errors  create  nonlinear 
coefficients  that  allow  noise  into  the  gradient  measurements. 


Figure  10:  Bell  Geospace  Air-FTG,  taken  from  [3]. 


The  Bell  Geospace  Air-FTG  is  a  3  disc,  rotating  accelerometer  type  gradient 
instrument  (shown  in  Figure  10)  that  was  launched  in  2002.  It  is  based  on  Lock- 
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heed  Martin’s  3D  FTG  and  includes  proprietary  post  mission  processing  upgrades 
for  improved  performance.  Each  disc  is  mounted  such  that  their  axes  of  rotation 
are  mutually  perpendicular  and  each  make  the  same  angle  with  the  vertical  (Fig¬ 
ure  11).  This  is  known  as  an  umbrella  conhguration.  The  GGIs  are  also  mounted  on 
a  three-gimballed  stabilizing  platform. 


Figure  11:  Bell  Geospace  Air-FTG,  taken  from  [3]. 

To  minimize  bias  from  the  orientation  or  movement  of  the  instrument,  the  as¬ 
sembly  of  rotating  discs  is  rotated  at  a  constant  rate  (300  deg/hr)  about  a  vertical 
axis.  The  Air-FTG  is  widely  used  in  airborne  gravity  gradient  mapping  in  the  USA, 
Ganada,  South  Africa,  Botswana,  Zambia,  and  Mali.  With  a  weight  of  roughly  A'bQkg 
and  requiring  approximately  1  cubic  meter  of  space  (with  data  acquisition  equip¬ 
ment),  the  Air-FTG  is  flown  primarily  in  the  Gessna  Grand  Garavan  -  though  it  has 
been  carried  by  zeppelins  (airships)  for  improved  stability  and  reduced  noise  [20] .  The 
Gessna’s  propeller  speeds,  engine  noise,  vibrations  and  other  disturbances  acting  on 
or  within  the  Air-FTG  are  monitored  during  each  flight  and  compensated  for  during 
post-flight  data  processing.  Since,  in  the  real  world,  no  instrument  is  perfect,  there 
is  some  non-linear  behavior  within  the  gradiometer.  These  nonlinear  coefficients  can 
cause  noise  due  to  host  vehicle  accelerations  and  disk  bearing  noise  within  the  de¬ 
sired  bandwidth.  This  noise  is  not  a  direct  measurement  of  host  vehicle  accelerations 
but  instead  is  the  various  products  of  acceleration  and  the  accelerometer  nonlinear 
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coefficients.  Thus,  if  the  coefficients  are  known  and  the  host  vehicle  accelerations  and 
measured,  the  noise  can  be  determined  and  eliminated.  This  is  exactly  what  BellGeo’s 
proprietary  High  Rate  Post  Mission  Compensation  (HRPMC)  does.  By  recording  the 
host  vehicle  accelerations  and  then  multiplying  these  accelerations  of  the  correct  order 
to  the  assumed  coefficients,  the  nonlinear  coefficients  are  found.  A  numerical  regres¬ 
sion  technique  is  then  performed  on  each  coefficient  until  the  noise  is  minimized. 
This  HRPMC  technique  has  been  proven  adequate  for  removing  noise  for  host  vehicle 
accelerations  of  around  O.lg  standard  deviation.  Note  that  two  other  factors  that 
can  induce  measurement  noise  are  misalignment  of  the  combination  of  accelerometers 
within  each  CCI  with  respect  to  the  plane  of  rotation  and  any  scale  factor  difference 
between  the  two  accelerometer  pairs.  Both  of  these  issues  are  addressed  before  each 
survey  through  calibration  techniques.  Additionally,  gravity  gradient  measurements 
are  extremely  sensitive  to  gravitational  held  disruptions  caused  by  nearby  masses. 
Such  masses  include  the  host  vehicle  structure  and  stores.  Since  these  masses  move 
with  the  host  vehicle,  it  is  critical  to  remove  their  inhuence  from  the  measured  data. 
This  process  is  accomplished  by  hying  a  specially  designed  survey  pattern.  Airborne 
gradiometer  surveys  are  always  conducted  in  an  orthogonal  pattern  which  results  in 
many  crossing  points.  Data  from  these  points  is  then  used  to  remove  host  vehicle 
gravitational  ehects  in  part  of  a  Low  Rate  Post  Mission  Compensation  (LRPMC) 
process  [21]. 

For  relatively  good  resolution,  surveys  are  typically  hown  using  drape  methods 
at  50-100m  above  the  ground  since  the  signal  strength  in  the  Air-FTC  drops  oh 
with  the  cube  of  the  distance  to  the  target.  The  current  resolution  of  the  Air-FTC, 
after  HRPMC  and  pre-hight  calibration,  is  approximately  5  Eotvos  at  a  gradient 
production  rate  of  IHz  with  a  spatial  resolution  of  several  hundred  meters.  Without 
the  aforementioned  processing  steps,  the  Air-FTC  noise  levels  are  approximately  12- 
IbEo  [3,21].  Note  that  raw  accelerometer  data  is  sampled  at  over  IQQHz  but  after 
compensation  and  demodulation  processes,  gradients  are  produced  at  IHz  [22].  It  is 
also  noted  that  gradiometer  noise  specihcations  are  often  given  in  terms  of  a  noise 
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spectral  density  (NSD)  having  nnits  of  Eoj ^ Hz.  Fnrthermore,  gravity  gradiometer 
manufactnrers  claim  zero-mean  Gaussian  white  noise  for  their  GGIs.  These  white 
noise  characteristics  are  only  valid  over  a  certain  range  of  frequencies  (or  bandwidth) 
-  generally  from  0  to  the  Nyquist  rate  or  a  low  pass  hlter  cutoff  frequency  [23]. 

The  Falcon  AGG  is  another  rotating  accelerometer  type  gravity  gradiometer. 
The  technology  was  jointly  developed  by  BHP  Billiton  and  Lockheed  Martin  and 
recently  sold  to  Fugro  NV.  The  Falcon  was  considered  the  hrst  airborne  gravity  gra¬ 
diometer  -  initially  flying  in  1997  and  used  for  survey  work  in  1999.  Fugro  has  suc¬ 
cessfully  used  the  Falcon  for  airborne  gravity  surveys  in  4  aircraft  (3  Gessna  Grand 
Garavans  and  1  helicopter).  This  system  measures  only  two  components  of  the  grav¬ 
ity  tensor  {T^z  and  Tyz)  and  uses  these  to  calculate  the  vertical  component  of  the 
tensor,  Tzz-  The  vertical  gravity  gradient  RMS  noise  is  around  5Eo  after  post  flight 
processing  techniques  similar  to  those  of  the  Bell  Geospace  Air-FTG.  A  6th  order 
Butterworth  hlter  with  a  cutoff  wavelength  of  400m  is  typically  used  in  Falcon  AGG 
data  processing  [24,25].  Assuming  a  survey  speed  of  60m/s,  this  corresponds  to  a 
cutoff  frequency  of  approximately  0.15Hz.  Dimensions  and  weight  of  the  Falcon  and 
its  data  acquisition  equipment  are  similar  to  those  of  the  Air-FTG. 

The  ARKeX  FTGeX  is  a  GGI  very  similar  to  the  previously  mentioned  instru¬ 
ments.  It  too  is  based  on  Lockheed-Martin’s  3D-FTG  and  is  often  used  with  ARKeX 
proprietary  technology  known  as  BlueQube.  BlueQube  involves  the  combination  of 
gravity  gradiometry,  magnetic  gradiometry,  digital  terrain  mapping  (LiDAR),  and 
digital  video  to  construct  a  complete  map  of  the  surveyed  area.  As  with  the  Air-FTG 
and  Falcon  AGG,  RMS  noise  of  the  FTGeX  is  around  5Eo  after  post-hight  processing, 
while  its  size  and  weight  are  also  similar  [26]. 

While  several  versions  of  the  rotating  accelerometer  GGIs  have  been  presented, 
one  hnal  point  about  this  type  of  gravity  gradiometer  must  be  made.  These  are 
the  only  type  of  gradiometers  successfully  used  in  airborne  surveys.  All  other  types 
discussed  herein  are  either  in  early  flight  test  or  a  laboratory  setting. 
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Superconducting  Gravity  Gradiometer.  Superconducting  Gravity  Gradiome- 
ters  (SGGs)  get  their  name  from  the  type  of  accelerometer  that  is  used  in  the  in¬ 
strument.  These  accelerometers  rely  on  the  Meissner  effect  and  flux  quantization  to 
levitate  a  proof  mass  and  measure  the  force  required  to  hold  that  mass  in  place.  Pairs 
of  superconducting  accelerometers  provide  gradient  measurements  with  low  noise  and 
high  resolution.  They  do  so  because  superconductivity  and  extremely  low  tempera¬ 
tures  naturally  give  low  noise,  negligible  scale  factor  drift  and  mechanical  stability. 
Superconducting  circuits  can  also  be  balanced  such  that  their  responses  to  gravity 
gradients  are  largely  independent  of  all  linear  and  angular  accelerations  applied  to 
the  instrument.  This  balance  stems  from  the  ability  to  regulate  currents  in  the  vari¬ 
ous  superconducting  loops.  It  is  because  of  this  balance  that  the  scale  factors  remain 
perfectly  stable  in  time  [27]. 

Diving  deeper  into  the  physics  behind  the  SGG,  the  accelerometer  itself  is  exam¬ 
ined.  First,  a  brief  overview  of  superconductivity  and  the  Meissner  effect  is  presented. 
A  superconductor  is  dehned  as  an  “element,  inter- metallic  alloy,  or  compound  that 
will  conduct  electricity  without  resistance  below  a  certain  temperature”  [28].  Meissner 
discovered  that  when  a  superconductor  is  placed  in  a  weak  external  magnetic  held, 
the  held  only  penetrates  the  superconductor  for  a  very  short  distance,  after  which 
it  drops  rapidly  to  zero.  In  essence,  a  superconductor  will  expel  all  magnetic  helds 
(time  variant  and  invariant).  A  magnet  moving  by  a  “normal”  conductor  induces  cur¬ 
rents  in  the  conductor.  This  is  the  principle  on  which  an  electric  generator  operates. 
But,  in  a  superconductor,  the  induced  currents  exactly  mirror  the  held  that  would 
have  otherwise  penetrated  the  superconducting  material  -  causing  the  magnet  to  be 
repulsed.  The  Meissner  ehect  is  so  strong  that  a  magnet  can  actually  be  levitated 
over  a  superconductive  material,  as  shown  in  Figure  12  [29]. 

Flux  quantization  is  a  quantum  phenomenon  in  which  the  magnetic  held  is 
quantized.  This  occurs  in  type  II  superconductors  subjected  to  a  magnetic  held. 
Type  II  superconductors  are  characterized  by  their  gradual  transition  from  the  su¬ 
perconducting  to  the  normal  state  as  temperature  increases.  They  tend  to  be  made 


Figure  12:  A  Magnet  Levitating  above  a  Superconductor  due  to  the  Meissner  Effect, 
taken  from  [30]. 

of  metal  alloys  or  complex  oxide  ceramics.  Below  a  temperature  dependent  critical 
magnetic  field  Hd,  all  magnetic  flux  is  expelled  according  to  the  Meissner  effect  and 
perfect  diamagnetism  is  observed.  Up  to  a  temperature  dependent  second  critical  field 
value,  Hc2,  flux  penetrates  in  discrete  units  while  the  bulk  of  the  material  remains 
superconducting.  Within  this  group  of  type  II  superconductors  are  high  tempera¬ 
ture  superconductors.  High-temperature  superconductivity  allows  some  materials  to 
support  superconductivity  at  temperatures  above  the  boiling  point  of  liquid  nitrogen 
(approx.  77°  Kelvin).  As  such,  they  offer  the  highest  transition  temperatures  of  all 
superconductors.  The  ability  to  use  relatively  inexpensive  and  easily  handled  liquid 
nitrogen  as  a  coolant  has  increased  the  range  of  practical  applications  of  superconduc¬ 
tivity  [31].  Unfortunately,  this  higher  temperature  of  operation  will  make  a  system 
intrinsically  noisier,  thus  high  temperature  superconductors  are  not  suited  for  gravity 
gradiometer  use  [27].  Instead,  SGGs  must  use  type  II  low  temperature  superconduct¬ 
ing  material  maintained  at  approximately  4°  Kelvin  in  the  circuits  (loops)  within  the 
accelerometers. 

With  superconductivity,  the  Meissner  effect,  and  flux  quantization  now  defined, 
an  example  of  a  superconducting  accelerometer  is  presented.  Suppose  a  time  varying 
current  is  passed  through  a  coil  outside  of  a  superconductor.  This  coil  will  send 
out  a  field  that  will  induce  a  surface  current  on  the  superconductor.  Noting  that 
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the  superconductor  will  expel  all  magnetic  fields  up  to  a  material  and  temperature 
dependent  critical  field  (i.e.  the  Meissner  effect),  this  surface  current  will  interact 
with  the  current  in  the  coil  to  produce  a  repulsive  force  between  the  two  objects. 
Using  image  theory,  it  can  be  shown  that  the  surface  current  on  the  superconductor 
is  equivalent  to  having  an  image  of  the  coil  within  the  superconductor  itself.  This 
image  coil  is  exactly  the  same  distance  below  the  surface  as  the  real  coil  is  above  it 
as  shown  in  Figure  13. 


Figure  13:  Coil  and  Image  Coil  near  a  Superconductor,  taken  from  [32]. 


Now  the  field  confined  between  the  coil  and  superconductor  is  calculated:  B  = 
HoTiti,  where  Ho  is  the  permeability  of  the  material  (i.e.  how  susceptible  it  is  to  being 
magnetized),  n*  is  the  turns  per  meter  in  the  coil  and  i  is  the  current  in  the  coil. 
Recognizing  that  the  total  magnetic  energy  in  the  system  is  the  field  energy  per  unit 
volume  times  the  volume  of  the  space  between  the  coil  and  superconductor: 


Magnetic  Energy 


B^  1 

- Acd  =  -  Acdi^ 

2fio  ^ 


(27) 


where  Ac  is  the  area  of  the  coil  and  d  the  distance  between  the  superconductor  and 
the  coil.  Recall  that: 


Magnetic  Energy 


(28) 
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where  L  is  the  inductance.  Combining  Equations  27  and  28  gives: 

L  =  fiou'^Acd  +  0{(f)  (29) 

Equation  29  shows  that  the  inductance  of  the  coil  is  proportional  to  its  separation 
from  the  superconductor  surface.  Now,  if  the  proof  mass  is  made  of  superconducting 
material  and  is  introduced  in  the  vicinity  of  the  coil,  a  repulsive  force  is  present  as 
long  as  current  is  flowing.  This  force  is  given  by: 

+  (30) 

where  the  stiffness  of  this  magnetic  spring  is  determined  by  coil  non-linearities. 
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Figure  14:  Superconducting  Accelerometer  Schematic,  taken  from  [32]. 

Now  a  system  where  a  closed  superconducting  loop  levitates  a  superconducting 
proof  mass  (Figure  14)  is  constructed,  noting  that  as  the  superconducting  loop  passes 
below  its  transition  temperature,  the  magnetic  flux  in  the  superconducting  loop  re¬ 
mains  absolutely  stable  (by  flux  quantization)  and  has  no  noise  on  it!  Whenever  the 
proof  mass  moves  due  to  a  change  in  acceleration,  the  coil  inductance  changes  and 
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the  current  must  change  to  preserve  the  original  flux  in  the  loop.  This  new  current 
now  exerts  a  different  force  on  the  proof  mass  in  order  to  preserve  the  mass’s  location. 
This  change  in  current  is  detected  by  a  device  called  a  SQUID.  A  SQUID  is  a  Su¬ 
perconducting  Quantum  Interference  Device  and  is  the  most  sensitive  sensor  known 
to  science.  It  is  used  to  measure  the  changes  in  magnetic  fields  from  which  changes 
in  currents  can  be  determined.  While  the  superconducting  accelerometer  is  a  sensor 
within  a  sensor,  the  SQUID’s  resolution  of  around  6.21  x  10-21  Wb{Webers)/^/Hz 
is  very  accurate  [33].  To  put  this  in  perspective,  the  SQUID  can  sense  changes  in 
magnetic  fields  approximately  16  orders  of  magnitude  smaller  than  that  produced  by 
a  small  refrigerator  magnet. 


mi 


Figure  15:  Superconducting  Gradiometer  Schematic,  taken  from  [32]. 

Now  that  the  superconducting  accelerometer  has  been  presented,  the  supercon¬ 
ducting  gravity  gradiometer  is  easily  shown.  If  two  masses  and  two  loops  are  used 
(Figure  15),  differential  movements  in  the  proof  masses  can  be  detected.  If  both  proof 
masses  move  the  same  distance  in  the  same  direction  (downwards  for  example),  flux 
quantization  for  each  loop  requires  that  U  and  I2  increase.  Since  these  currents  flow 
in  opposite  directions  through  the  inductor  next  to  the  SQUID,  they  cancel  each  other 
out  and  the  SQUID  measures  no  change.  However,  if  the  masses  move  in  different 
directions  from  one  another  (or  if  one  moves  and  the  other  does  not),  Ji  and  I2  will 
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be  different  and  the  SQUID,  because  of  its  location,  will  pick  up  the  difference  in 
current.  The  current  difference  due  to  differential  movement  of  the  test  masses  is  a 
direct  result  of  a  gravity  gradient  [32], 

The  first  superconducting  GGI  was  developed  by  Dr.  Ho  Jung  Paik,  a  Univ. 
of  Maryland  professor,  in  conjunction  with  the  NASA  Goddard  Gryogenics  branch. 
This  hrst  version  of  the  SGG  consisted  of  three  orthogonal  pairs  of  superconducting 
accelerometers,  each  capable  of  measuring  linear  accelerations  and  the  gravity  gradient 
in  all  three  axes.  These  superconducting  accelerometers  are  made  of  high  purity 
niobium  (one  of  three  Type  II  superconducting  elements)  that  are  kept  cool  in  a 
Helium  bath.  Pancake  coils  are  used  to  levitate  the  proof  masses  which  are  suspended 
initially  by  a  weak  spring  [34]. 

An  equivalent  gravity  gradient  noise  of  Eotvos/ \/Hz  has  been  demon¬ 

strated  in  a  laboratory.  While  this  is  much  better  than  rotating  accelerometer  GGIs, 
the  Univ.  of  Maryland  superconducting  gravity  gradient  instrument  (SGGI)  fell  short 
in  intrinsic  noise  levels.  The  largest  contributors  to  these  noise  levels  were  the  earth’s 
gravitational  held  and  simulated  host  vehicle  acceleration  coupling  into  the  gradient 
outputs  through  various  mechanical  errors.  The  instrument  was  also  sensitive  to  ther¬ 
mal  huctuations  of  the  helium  bath,  liquid  helium  motion  and  boiloff,  and  particle 
heating.  Most  errors  could  be  controlled  by  precise  design  and  alignment  of  the  in¬ 
strument  or  removed  by  measuring  the  disturbances  (similar  to  techniques  used  in 
the  rotating  accelerometer  GGIs).  While  NASA  had  interest  in  testing  the  SGG  in 
space,  funding  was  not  available  for  such  a  mission  [34]. 

Updates  to  the  SGGI  have  involved  making  it  operable  in  a  moving  base  envi¬ 
ronment  and  include  revised  angular  accelerometers  for  measuring  the  gradient  along 
a  single  axis  while  using  three  linear  accelerometers  to  correct  for  residual  coupling 
to  linear  acceleration  due  to  imperfect  mass  balances.  This  device,  shown  in  Figure 
16,  is  known  as  the  UM-SAA  (University  of  Maryland  Superconducting  Angular  Ac¬ 
celerometer).  To  eliminate  the  thermal  sensitivity  issues,  the  updated  SAA  is  cooled 
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Figure  16:  University  of  Maryland  SAA,  taken  from  [35]. 

by  a  closed  cycle  refrigerator  based  on  a  dual-stage  pulse-tube  cold-head  known  as 
a  cryostat  (Figure  17).  The  pulse-tube  has  no  reciprocating  piston  in  the  cold-head 
thereby  greatly  reducing  the  harmonics  of  pressure  pulses  [35]. 


Figure  17:  University  of  Maryland  Cryostat,  taken  from  [35]. 

The  most  recent  University  of  Maryland  airborne  SAA  has  an  estimated  error  of 
approximately  0.3Eo  with  a  gradient  production  rate  of  IHz  [35].  While  bandwidth 
specihcations  are  not  clearly  stated  in  open  literature,  Lumley  et  ah  [32]  cite  “best 
spatial  resolutions  of  a  few  hundred  meters” . 

The  High  Density  Airborne  Gravity  Gradiometer  (HD-AGG)  is  a  another  air¬ 
borne  superconducting  GGI  currently  in  test  [36].  Greated  by  researchers  from  the 
University  of  Maryland,  University  of  Western  Australia,  Ganadian  Space  Agency 
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and  Gedex,  the  HD-AGG  is  a  gradiometer  designed  to  be  carried  in  typical  geolog¬ 
ical  snrvey  aircraft  (likely  the  Gessna  Grand  Garavan)  and  has  similar  dimensions 
to  airborne  rotating  accelerometer  GGIs.  Using  licensed  technology  from  the  Gana- 
dian  Space  Agency,  Gedex  added  an  external  distnrbance  isolation  platform,  known 
as  the  GeoMIM,  to  the  sensor.  Similar  isolation  platforms  are  fonnd  on  the  Space 
Shnttle  [37].  The  HD-AGG  has  reportedly  achieved  a  gravity  gradient  error  variance 
of  less  than  lEo  for  measurements  made  every  second  -even  in  moderate  turbulence 
(Im/s^).  Main  [38]  cites  spatial  resolutions  of  60m  at  hxed  wing  survey  aircraft 
speeds.  As  such,  De  Beers  entered  a  strategic  agreement  with  Gedex  in  2006  for  use 
of  the  HD-AGG  in  diamond  detection  [39]. 


Rounding  out  the  list  of  superconducting  GGIs  currently  in  test  is  the  ARKex 
Exploration  Gravity  Gradiometer  (EGG).  The  EGG,  shown  in  Figure  18,  was  de¬ 
veloped  by  Dr.  John  Lumley  at  Oxford  Instruments  Superconductivity  Ltd  with 
assistance  from  ARK  Geophysics  Ltd  (now  a  part  of  ARKeX)  and  was  set  to  enter 
a  trial  deployment  sometime  in  2008  [26].  Though  exact  details  of  the  EGG  are  not 
published,  it  likely  uses  technology  from  the  UM-SAA  and  the  University  of  Western 
Australia’s  Orthogonal  Quadrupole  Responder  (UWA  OQR).  Data  from  lab  testing 
indicates  similar  error  variance,  noise  properties,  and  gradient  production  rates 
as  the  HD-AGG  [12]. 


Figure  18:  ARKeX  EGG,  taken  from  [40]. 
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Atom  Interferometer  Gravity  Gradiometer.  Atom  Interferometer  GGIs  are 
based  on  quantum  particle-wave  duality.  According  to  quantum  mechanics,  atoms 
behave  like  waves,  as  does  light.  Therefore,  an  interferometer  (a  device  that  shows 
the  pattern  of  interference  created  by  the  superposition  of  two  or  more  waves)  that 
examines  the  properties  of  these  atom-waves  can  be  constructed.  Because  atoms  have 
a  hnite  mass,  they  are  extremely  sensitive  to  changes  in  gravity.  In  an  Atom  Interfer¬ 
ometer  accelerometer,  “beams”  of  atoms  are  split  into  two  parts  via  a  beamsplitter 
and  then  allowed  to  travel  a  hnite  distance.  If  gravity  is  acting  on  these  atoms  as  they 
travel  over  a  certain  length  (i.e.  the  atoms  are  under  the  inhuence  of  gravity),  the 
interferometer  will  pick  up  a  phase  shift  that  will  affect  the  phase  and/or  frequency 
of  the  measurements  [41]. 

The  Stanford  University  Atomic  Interferometer  Gradiometer  is  a  joint  effort 
between  Stanford  University  and  the  NASA  Jet  Propulsion  Laboratory.  It  is  cur¬ 
rently  a  ground  based  sensor  that  will  soon  be  tested  as  part  of  an  all-atom,  gravity- 
compensated  inertial  navigation  system.  The  sensor  will  be  mounted  inside  an  RV 
and  driven  around  the  country  to  measure  INS  drift  under  realistic  conditions.  The 
gradiometer  itself  uses  two  quantum  gravity  accelerometers  located  a  certain  hxed 
distance  apart  (Figure  19).  Inside  these  accelerometers,  cooled  Gesium  atoms  are 
condensed  into  a  small  cloud  in  a  magneto-optic  trap  (MOT).  The  MOT,  shown  in 
Figure  20,  consists  of  three  pairs  of  counter-propagating  laser  beams  along  three  axes 
centered  about  a  non-uniform  magnetic  held  and  can  collect  up  to  10®  atoms  [41]. 

After  these  atoms  are  collected,  further  cooling  slows  them  to  an  RMS  velocity 
of  a  few  cm/s.  The  cold  atoms  are  then  launched  vertically  into  an  atomic-fountain 
so  that  the  sensors  have  twice  the  available  interaction  time  with  the  atoms  for  a 
given  height.  The  atom  interferometer,  shown  in  Figure  21,  is  made  up  of  Raman 
transitions  between  two  hyperhne  ground  states  with  a  f  —  tJ"  —  f  pulse  sequence. 
The  hrst  pulse  creates  an  equal  superposition  of  atoms  in  two  hyperhne  ground  states 
(beam  splitting).  The  second  and  third  pulses  redirect  and  recombine  the  atom- wave. 
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Figure  19:  Schematic  of  the  Stanford/ JPL  AI  Gradiometer,  taken  from  [42], 


Atom  Cloud  in 
Optical  Molasses 

Figure  20:  Magneto-Optic  Trap 

If  gravity  is  acting  on  the  atoms,  their  paths  will  be  different  and  a  phase  shift  will 
occur. 

This  phase  shift  is  given  by  A0  =  2krgT^,  where  Tj  is  the  interrogation  time 
(the  time  between  light  pulses),  kr  is  the  Raman  laser  wave  number  and  g  is  gravity. 
Gravity  gradients  are  sensed  and  quantihed  when  there  is  a  mismatch  between  the 
readings  of  two  hxed  position  accelerometers.  Laboratory  testing  of  the  Stanford/JPL 
AI  GGI  has  shown  gravity  gradient  sensitivity  of  10Eo/\^H z,  though  improvements 
are  expected  [41].  It  should  also  be  noted  that  the  gradiometer  is  approximately 
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Figure  21:  Atom  Interferometer 

1.25m  tall  (Figure  22)  and  also  requires  a  fairly  bulky  array  of  control  and  laser 
frames,  though  efforts  are  being  made  to  compact  the  system. 


Figure  22:  Stanford/JPL  AI  Gradiometer  Dimensions 
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Navigation  and  Terrain  Avoidance  via  Gravity  Gradiometry  -  Previous 
Works 

The  use  of  GGIs  as  navigation  aids  has  been  investigated  since  the  1960s.  Initial 
research  efforts  focused  on  real-time  measurement  of  the  gravity  anomaly  to  provide 
improved  unaided  inertial  navigation  accuracy.  Up  until  recently,  there  has  been 
relatively  little  gravity  gradiometer  map-matching  based  aircraft  navigation  or  ter¬ 
rain  avoidance  research  published  in  open  literature.  In  the  mid  1970s,  Metzger  and 
Jircitano  [43]  presented  an  investigation  into  using  gravity  and  gravity  gradient  map¬ 
matching  to  update  an  INS.  Host  vehicle  velocities  of  up  to  approximately  250m/s 
were  examined  and  it  was  found  that  gravity  gradients  provided  better  results  due  to 
better  signal  uniqueness. 

In  the  mid  1980s,  Bell  Aerospace  Textron  began  development  on  a  system  for 
enhanced  passive  submarine  navigation.  The  system  was  to  also  provide  real-time 
underwater  terrain  maps  in  areas  where  accurate  terrain  and  obstacle  data  may  be 
limited.  Around  that  time,  Affleck  and  Jircitano  [44]  proposed  an  INS  that  received 
position  updates  from  a  full  tensor  gradiometer/map-matching  algorithm.  The  study 
was  carried  for  an  aircraft  with  a  velocity  of  lOOm/s  and  altitudes  that  ranged  from 
100-400m  with  promising  results.  A  ship  based  navigation  performance  analysis  was 
executed  as  well,  also  with  encouraging  results.  However,  very  little  detail  about  the 
navigation  algorithm  was  given  -  likely  due  to  the  proprietary  nature  of  the  work.  In 
1994,  Jircitano  and  Dosch  [45]  patented  a  Gravity  Aided  INS  (GAINS)  using  a  GGI 
and  a  vertical  gravimeter  for  covert  submarine  navigation.  Shortly  thereafter.  White 
and  Jircitano  patented  “gradiometer  based  terrain  estimation”  [46].  The  system(s) 
came  to  fruition  in  the  mid  1990s  and  is  known  as  the  Lockheed- Mart  in  Universal 
Gravity  Module  (UGM).  The  UGM  consists  of  gravimeters  and  gravity  gradio meters 
and  implements  a  gravity  map- matching  algorithm  to  passively  bound  INS  errors.  It 
can  also  provide  real-time  underwater  terrain  maps  based  on  estimation  techniques 
applied  to  existing  databases.  In  1998,  the  UGM  was  successfully  tested  on  the  USS 
Memphis  fast  attack  submarine  [47].  Likely  due  to  the  covert  nature  of  the  business. 
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more  recent  information  on  the  UGM  or  Jircitano’s  work  was  unable  to  be  found  in 
open  literature. 

Archibald  [48]  took  a  different  map-matching  approach  for  his  PhD  disserta¬ 
tion  by  using  a  neural  network  based  map-matching  algorithm  to  correlate  gravity 
gradiometer  measurements  to  an  existing  map.  While  his  gradiometer  model  was  rel¬ 
atively  simple,  his  method  matched  points  among  large  amounts  of  geophysical  data 
and  could  be  used  as  an  acquire  mode  in  a  staged  map-matching  scheme. 

In  1995,  Gleason  [23]  presented  a  method  to  optimally  generate  gravity  gradient 
maps  and  discussed  the  effects  of  gradiometer  hltering  in  a  terrain  avoidance  scenario 
as  well  as  many  other  practical  issues  of  a  GGI-based  map  matching  scheme.  His 
work  was  the  hrst  found  to  address  GGI  sampling  rate,  gradient  production  rate, 
noise,  and  bandwidth  in  detail. 

Blaylock  et  ah  [49]  then  presented  a  terrain  avoidance  method  using  a  gravity 
gradiometer  (theoretically)  on  board  an  F-16.  Likely  an  extension  of  Jircitano’s  work, 
the  paper  did  not  go  into  detail  about  an  actual  terrain  avoidance  algorithm.  It  did, 
however,  present  some  estimated  gradiometer  performance  requirements  and  various 
GGI  signal  levels  as  modeled  obstacles  were  approached.  Additionally,  very  little 
information  was  given  on  solving  the  inverse  problem  of  calculating  obstacle  range 
from  gravity  gradients. 

Though  he  did  not  use  a  map-matching  algorithm,  Jekeli  [50]  showed  that  fu¬ 
ture  high  accuracy  IMUs  could  provide  near  GPS  accuracy  (5m  error  after  1  hour 
of  dead  reckoning)  if  the  gravity  error  was  compensated  with  a  full  tensor  gravity 
gradiometer  providing  Is  updates  with  O.li^o  of  RMS  noise.  The  premise  is  that 
if  highly  accurate  accelerometers  and  gyroscopes  are  used,  errors  due  to  bias,  scale 
factor/misalignments,  platform  tilt,  white  noise  and  random  walk  become  very  small 
relative  to  errors  from  uncompensated  gravity.  His  research  is  an  integral  part  of 
the  DARPA  Precision  Inertial  Navigation  Systems  (PINS)  program.  While  this  is 
arguably  the  more  elegant  and  simpler  approach  to  GGI  based  navigation,  dead  reck- 
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oiling  with  this  system  still  produced  levels  of  position  error  that,  over  time,  could 
hamper  mission  effectiveness.  However,  Jekeli  also  explains  the  drawbacks  of  matching 
GGI  measurements  to  a  pre-existing  map.  He  states  that  worldwide  map  coverage, 
particularly  in  remote  or  mountainous  areas  is  very  limited.  Also,  the  accuracy  of 
gravity  gradient  maps  is  derived  from  major  surveying  efforts  that  are  not  easily  aug¬ 
mented  when  improved  accuracy  is  required.  Finally,  he  cites  the  largest  map-making 
hurdle  to  overcome:  “the  required  gravitational  accelerations  are  the  horizontal  com¬ 
ponents  of  the  gravitational  vector  at  altitude  (for  aircraft  navigation),  whereas  the 
data  typically  are  vertical  components  on  the  Earth’s  surface  (being  the  most  easily 
measured)”  [50].  The  concerns  are  similar  to  those  expressed  in  the  early  days  of  TRN 
(such  as  TERGOM)  [2].  If  the  map  making  issues  mentioned  herein  are  solved,  the 
following  question  is  posed:  If  a  gradiometer  is  providing  INS  gravity  compensation 
and  is  already  present  within  the  system,  why  not  use  it  to  give  map-matching  based 
updates  to  the  INS  as  well?  This  could  potentially  give  the  best  of  both  worlds  -  low 
INS  error  under  dead  reckoning  and  the  ability  to  correct  the  growing  error  over  time. 

Most  recently,  Richeson  [51]  provided  an  in-depth  study  of  passive  navigation  via 
gravity  gradient  map-  matching  by  developing  an  INS  model  that  used  an  Extended 
Kalman  Filter  (EKF)  to  integrate  gravity  gradiometer  and  gradient  map  information 
in  a  position-updating  algorithm.  He  showed  that  a  gradiometer  providing  Is  updates 
with  O.OOlEo  of  noise  allowed  a  map-matching  algorithm  to  meet  GPS  performance 
levels.  While  his  research  neglected  terrain  effects  by  focusing  mainly  on  high  altitude 
(20A;m)  and  high  velocity  {2km/ s)  scenarios  seen  by  a  hypersonic  vehicle,  an  estima¬ 
tion  of  several  important  GGI  signal  threshold  requirements,  presented  in  Ghapter  3, 
were  deduced  from  his  work.  A  schematic  of  Richeson’s  navigation  system  is  shown 
in  Figure  23. 

This  research  aims  to  advance  knowledge  and  understanding  of  GGI-based  air¬ 
craft  navigation  and  terrain  avoidance  by  more  rigorously  modeling  the  GGI  and 
examining  the  signal  itself  over  a  wider  range  of  flight  conditions  than  previously 
studied.  Note  that  the  navigation  portion  of  this  study  will  focus  on  signal  usefulness 
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Figure  23:  Richeson’s  Proposed  Navigation  System,  taken  from  [51]. 


in  map-matching  scenarios  rather  than  utility  for  INS  gravity  error  compensation. 
Additionally,  terrain  avoidance  scenarios  involving  a  variety  of  hazardous  obstacles 
will  be  examined  in  search  of  a  useable  impact  warning  threshold. 
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III.  Methodology 


Overview 

To  determine  the  feasibility  of  GGI  based  passive  navigation  and  terrain  avoid¬ 
ance,  a  buildup  approach  to  simulating  the  GGI  signal  will  be  accomplished,  followed 
by  characterization  of  the  signal  level  needed  to  successfully  navigate  or  avoid  terrain, 
and  ultimately  culminating  in  a  comparison  between  required  and  achieved  signal 
levels.  First,  it  is  desired  to  bound  the  performance  metrics  of  GGI  based  navigation 
and  terrain  avoidance  by  simulating  the  GGI  signal  in  the  best  and  worst  case  sce¬ 
narios.  In  a  basic  sense,  these  scenarios  will  be  addressed  by  using  areas  with  rapid 
terrain  changes  and  areas  with  relatively  low  levels  of  terrain  changes  (i.e.  mountain¬ 
ous  versus  flat  terrain).  A  realistic  representation  of  the  earth’s  gravity  gradients  in 
these  test  areas  will  then  be  found  via  a  combination  of  modeling  techniques.  Once 
acceptable  maps  have  been  generated  for  all  test  conditions,  the  simulated  GGI  signal 
will  be  determined  by  manipulating  values  from  these  maps  with  appropriate  hlter- 
ing  and  noise  which  mirror  current  or  near-future  GGI  performance.  Finally,  the 
GGI  signal,  it’s  ability  to  be  correlated  to  the  original  gradient  maps,  and  ultimately, 
it’s  usefulness  for  navigation  will  be  investigated  qualitatively  and  by  metrics  based 
on  previous  works.  A  simple  case  study  on  the  threshold  signal  levels  indicating  an 
imminent  terrain  collision  will  be  accomplished  as  well.  All  work  herein  will  be  ac¬ 
complished  in  Matlab®  using  the  Microsoft  Windows  Vista  64  operating  system  on  an 
Intel  Gore2Duo  processor  at  3GHz  with  8GB  of  RAM.  All  code  used  in  this  research 
can  be  found  in  Appendix  A. 

Gravity  Gradient  Maps 

Before  the  signal  produced  by  the  GGIs  can  be  simulated,  representative  maps 
of  gravitational  gradients  produced  by  the  earth  must  be  computed.  The  method 
chosen  for  map  generation  is  the  combination  of  gradients  derived  from  the  Earth 
Gravitational  Model  1996  (EGM96)  and  gradients  derived  from  a  frequency  domain 
based  technique  similar  to  the  rectangular  prism  method  shown  in  Ghapter  1.  The 
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EGM96  gradients  will  be  used  to  account  for  long  wavelength  gravitational  effects 
that  generally  correspond  to  anomalies  well  beneath  the  earth’s  visible  topography 
while  the  rectangular  prism-like  method  will  be  used  to  model  terrain  effects.  The 
total  gravitational  gradient  disturbance  model  used  in  this  study  is  given  by  Tij  = 
T- 

’!'JEGM96  '  ’>'3  Terrain' 

It  must  be  stressed  that  these  methods  do  not  give  perfect  values  for  gravity 
gradients  (nor  do  any  methods  at  this  time)  but  can  be  calculated  relatively  quickly 
and  do  represent  realistic  trends  that  should  be  seen  by  a  GGI.  It  should  also  be 
noted  that  Tij  (the  disturbing  potential)  and  Uij  (the  normal  potential)  can  be  dehned 
differently  depending  on  the  application.  For  example,  the  oil  and  mining  industries 
will  often  generate  a  terrain  gradient  model  based  on  the  assumption  of  constant 
terrain  density  and  treat  it  as  a  part  of  Uij  in  order  to  subtract  it  out  and  find 
density  changes  within  the  terrain.  For  navigation  purposes,  terrain  generated  gravity 
gradient  information  is  a  requirement.  Knowledge  of  what  corrections  have  been 
applied  to  map  should  be  obtained  when  using  gradient  data  from  an  outside  source. 
For  example,  if  a  GGI  based  map-matching  navigation  system  database  were  loaded 
with  gravity  gradient  data  that,  unbeknownst  to  the  user,  had  terrain  effects  removed, 
the  system  would  be  rendered  useless. 

Earth  Gravitational  Model  1996.  As  previously  mentioned,  gravitational 
potential  outside  of  the  attracting  masses  follows  Laplace’s  equation  [6].  Spherical 
harmonics  are  the  angular  portion  of  an  orthogonal  set  of  solutions  to  Laplace’s  equa¬ 
tion  (represented  in  spherical  coordinates).  This  set  of  solutions  is  linear,  thus  the 
gravitational  potential  (or  disturbing  potential)  may  be  modeled  as  some  truncated 
series  of  spherical  harmonics  given  below  [52]: 


T{r,(l),X)  =  -j-  ^  ^  {Cnm  COSm\  +  Snm  siu  TuX)  P nm{cOS  (p)  (31) 


ae  ^ '  \r  J 

n=2  m=0 
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where  r,  0,  A  are  the  geocentric  distance,  lattitude  and  logitude,  respectively,  GM  is 
the  prodnct  of  the  gravitational  constant  and  the  earth’s  mass,  ae  is  the  semi-major 
axis  of  the  reference  ellipsoid,  N^ax  is  the  maximnm  degree  of  the  spherical  harmonic 
expansion,  n,m  are  the  degree  and  order,  Cnmt  Snm  are  the  normalized  geopotential 
coefficients  and  Pnm{cos6)  is  the  normalized  associated  Legendre  fnnction  (Legendre 
fnnctions  are  canonical  solntions  to  the  general  Legendre  differential  equation  that  is 
encountered  when  solving  Laplace’s  equation  in  spherical  coordinates). 

To  date,  there  are  a  variety  of  global  geopotential  models  which  express  the 
Earth’s  potential  field  in  terms  of  spherical  harmonic  coefficients.  These  models  are 
derived  from  satellite  orbit  tracking,  terrestrial  gravimetry,  satellite  altimetry,  or  air¬ 
borne  gravimetry  (or  a  combination  of  these  and  other  methods)  and  used  to  compute 
a  gravimetric  geoid  [53] .  The  geoid  (as  it  pertains  to  earth)  is  defined  as  the  equipo- 
tential  surface  of  the  earth’s  gravity  held  coinciding  with  the  mean  sea  level  (MSL)  of 
the  oceans  [54].  In  very  broad  terms,  the  geoid  is  a  mathematical  hgure  of  the  earth’s 
surface  dehned  by  gravitational  measurements,  as  opposed  to  the  smooth  surface  of 
a  reference  ellipsoid  such  as  the  WGS84  ellipsoid  (see  Figure  24).  In  other  words, 
it  is  a  surface  which  best  hts  the  mean  sea  level  without  winds,  ocean  currents,  or 
other  disturbing  forces.  While  the  geoid  and  ellipsoid  surfaces  end  up  being  similar  in 
practice,  the  geoid  surface  varies,  or  undulates,  approximately  -|-85m  to  -106m  with 
respect  to  the  WGS84  reference  ellipsoid  and  can  change  slightly  depending  on  the 
method  used  to  calculate  the  geoid. 


Topographic 

Surface 


Figure  24:  Exaggerated  Illustration  of  the  Geoid,  Ellipsoid,  and  Topography. 
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One  of  the  most  commonly  used  models  is  EGM96.  The  EGM96  geopotential 
model  is  a  composite  solution  consisting  of  spherical  harmonic  coefficients  to  degree 
and  order  360  (n  and  m  in  Equation  31,  respectively).  EGM96  is  made  up  of  data  from 
various  contributors  and  completed  by  a  joint  effort  between  the  National  Imagery  and 
Mapping  Agency,  NIMA  (now  the  National  Geospatial-Intelligence  Agency,  NGA), 
the  NASA  Goddard  Space  Flight  Genter  and  The  Ohio  State  University.  Some  of 
the  data  sources  include:  gravity  data  from  NIMA  obtained  by  airborne  surveys  and 
other  gravity  collection  processes,  data  from  the  GEOSAT  Geodetic  Mission  (a  US 
Navy  satellite  with  a  RADAR  altimeter  capable  of  measuring  distances  to  the  sea 
surfaces  within  5cm),  and  data  from  the  ERS-1  satellite  [55]. 

The  resolution  of  a  particular  model  is  given  by  'kR/u  where  R  is  the  earth’s 
average  radius  and  n  is  the  harmonic  degree  of  the  model.  Using  R  =  6371000  m,  the 
spatial  resolution  of  the  EGM96  model  is  approximately  56  km.  Since  the  resolution 
of  the  EGM96  and  rectangular  prism  maps  often  differ,  a  spline  interpolation  is  often 
applied  when  htting  the  EGM96  data  to  other  grids.  To  better  illustrate  the  shape  of 
this  particular  geoid,  EGM96  undulations  for  a  tide-free  system,  with  respect  to  the 
WGS84  ellipsoid,  are  shown  in  Figure  25. 
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Figure  25:  EGM96  Geoid  Undulations  with  respect  to  the  WGS84  Ellipsoid. 
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The  five  independent  EGM96  gravitational  disturbance  gradients  were  calcu¬ 
lated,  using  the  parameters  listed  in  Table  1,  via  modified  freeware  using  the  following 
relationships  [53,56]: 


XX  EGM9Q  ^^3 


max 

^  ^  ^  ^  ^  ^  J  (^nm  COS  Tfl\  -\-  Sum  StTlTTlX^  X  (^Clj2rn,-^n,m— 2 

n=2  m=0 


T. 


“1“  [^nm  “1“  S'  “h  ^72771-^72,771+2  )  7 

(32) 

(S ^  (XC  77  I  3 

^  ^  ^  ^  ^  77777  COS  TOX  S^im  SiflTTlX^  X  (c?72r77-P77— 1,771  —  2 


77=2  777=1 


T, 


“1“  5^77777-^77— l,777  (^^^0)  “1“  ^77777-^77— 1,777+2  )  5 

(33) 

(S  ]}d  ^  QjC  77  I  3 

^  ^  ^  ^  77771  COS  TTlX  -\-  Sfim  sin  TTl/\)  X  (/3j2777-P77,771— 1 


-^^GM96  ^^3  V  r  7 

77=2  771=0 

“1“  V^77771-P77,777+l  )  5 

(34) 

^  /oe\»+3,^  ,  ,  . 

Xyy  EGM9Q  3~  /  ^  ^  (  1  (^77771  COS  TflX  “h  bnm  StnTflXj  X  (^(Xj^tti-^  77,771— 2  (^^^0) 

77=2  771=0 

H“  bnmPnm{'^i‘Cf^4^)  S'  0^2772 -Pri, 7n+2 ('^XTXfJ^) ) , 

(35) 

^  ^  ^  ^  {G nm  COSTTlX  Snm  siflTTlX^  X  (sf?r0) 


T, 


2/^i.GM96  A^\r  J 

n=2  m=l 


T  hnm-Pn— l,m-|-l('S^'^0))) 


(36) 
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where  a,  b,  c,  d,  g,  h,  (3,  /i  and  g  are  the  coefficients  of  the  Legendre  functions 
and  are  located  in  Appendix  B.  Table  1  summarizes  the  properties  used  in  the  EGM96 
based  gradient  calculations. 


Table  1:  EGM96  Parameters 


Parameter 

Value 

GM 

3986004.415  x  lO^m^s^ 

a 

6378136.3  m 

G2,0 

Tide  Free 

Reference  Ellipsoid 

WGS84 

It  should  be  noted  that  the  newer  EGM2008  model  was  considered  for  computa¬ 
tion  of  the  gradients  due  to  its  degree  and  order  of  2159  (and  corresponding  resolution 
of  approximately  9km)  [57].  However,  due  to  the  recursive  generation  used  for  the 
associated  Legendre  functions,  these  functions  can  become  unstable  at  higher  degrees 
(approximately  2100).  Even  with  algorithms  that  allow  for  better  stability  at  high 
degrees,  preliminary  results  have  shown  that  the  estimation  of  gravity  components 
can  take  a  considerable  amount  of  time  [53].  Within  the  scope  of  this  investigation, 
the  drawbacks  of  using  the  EGM2008  model  outweighed  the  increased  resolution. 


Extended  Parker  Method.  Since  terrain  effects  make  up  the  largest,  most 
rapidly  changing  part  of  the  bias  removed  GGI  signal  at  lower  altitudes  [23],  an 
approximation  of  these  effects  must  be  modeled  to  get  a  relatively  accurate  signal 
representation  for  many  of  the  simulation  test  points.  Nagy’s  formulae  (Equations 
8-12)  for  determining  gradients  at  a  point  from  a  single  rectangular  prism  can  be 
expanded  to  include  the  effects  of  an  entire  grid  of  rectangular  prisms  as  follows: 


T, 


NM 


total 


n=\ 


(37) 


where  NM  corresponds  to  the  number  of  row  and  column  entries  on  the  grid.  While 
proper  use  of  the  rectangular  prism  method  gives  good  insight  into  the  general  behav¬ 
ior  of  the  terrain  implied  gradients,  it  can  be  extremely  expensive  computationally 
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for  a  large  grid.  Parker  [58]  presents  a  fast  frequency  domain  method  of  calculating 
the  potential  while  Jekeli  and  Zhu  [59]  apply  it  to  the  computation  of  gravitational 
gradients. 

If  the  surface  of  the  geoid  is  approximated  as  a  plane  and  a  constant  density 
contrast  is  assumed,  the  gravitation  potential  due  to  the  volume  mass  (terrain)  of 
height  h,  bounded  beneath  by  the  area  A  (approximated  as  a  Om  MSL  plane),  can  be 
written  in  a  form  similar  to  Equation  1: 

h 

V  =  GAp  jjj\  dz'dA  (38) 

A  0 


where  h  =  h{x',y')  and  is  the  terrain  height  at  each  point. 

According  to  2-D  Fourier  transform  theory,  if  g{x,y)  is  a  hnite  energy  function: 


{g{x,y)Ydxdy  <  oo 


(39) 


—  OO  — oo 


Then  there  exists  a  2D  continuous  Fourier  transform  pair: 


G{h,f2)  =  ^{g{x,y))=  j  y  (7(a;,|/)xe-*2-(h-+A^)da;d|/ 

—  OO  — OO 

oo  oo 

g(a:,9)  =  !a-‘(G(/i,/2))=  f  /  G(/i, /J  X 


—  OO  — oo 


(40) 

(41) 


where  /i  and  /2  are  spatial  frequencies  corresponding  to  coordinates  x  and  y  and  A  de¬ 
notes  a  2-D  Fourier  transform.  Physically  speaking,  the  Fourier  transform  is  a  method 
to  break  a  function  into  oscillatory  components  (i.e.  a  frequency  domain  representa¬ 
tion).  The  inverse  Fourier  transform  sends  the  function  from  its  frequency  domain 
representation  back  into  the  spatial  domain.  This  technique  is  particularly  benehcial 
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because  some  computations  (such  as  convolution)  can  be  accomplished  much  faster 
in  the  frequency  domain. 

If  the  assumption  is  made  that  A  spans  from  —  cx)  to  cx)  in  two  dimensions  is 
made,  the  Fourier  transform  of  V  becomes  [59] : 

h 

^  dz'dA  (42) 

A  0 


A(I/)  =  GAp 


which,  after  using  polar  coordinates,  can  be  expressed  as: 


-2-Kfz 

2nP 


('g27r/h(x',y')  _  ]^^g-i27r(/ix'+/2F)^^ 


(43) 


where  /  =  a/ /f  +  /| .  By  expanding  in  Equation  43  via  a  Taylor  series, 

A(E)  becomes: 

OO 

3(1/)  =  27TGe\pe-^’l’  Y1  -(2’//)"^"  y'))")  (44) 

1 

n=l 


To  obtain  the  gradients  from  the  potential,  a  frequency  domain  relationship  provided 
by  Jekeli  [60]  is  used: 


A(Id,)  =  /i,,A(E) 


(45) 


Recall  that  V  was  dehned  as  the  potential  due  to  some  arbitrary  mass  volume.  If  the 
mass  volume  is  defined  as  the  terrain  above  the  geoid.  Equation  45  becomes: 

terrain)  ~  Pij‘A{T terrain)  (46) 
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where: 


/ill  =  -(27^)Vl^  hi2  =  -(27r)Vi/2,  /ii3  =  -i(27r)Vi/  , 

/i22  =  -(27r)V|  ,  h23  =  -*(27r)V2/  , 

h33  =  (27r)V^  (47) 

By  applying  an  inverse  Fourier  transform  to  the  frequency  domain  representation  of 
Tij  in  Equation  46,  the  terrain  implied  gravity  disturbance  gradients  are  found: 

r., 2irGAp3-‘  (48) 

In  practice,  some  assumptions  regarding  Equation  48  must  be  made.  First,  A  is  a 
hnite  area  corresponding  to  the  area  of  the  elevation  grid  thus  Equation  48  becomes  an 
approximation.  Also,  the  assumption  is  made  that  Fourier  transforms  of  the  powers 
of  h{x',  y')  exist.  Since  h  is  given  at  discrete  points  in  the  elevation  grid  (Figure  A.l), 
discrete  approximations  of  the  continuous  Fourier  transform  are  used.  Furthermore, 
a  hnite  Taylor  series  expansion  is  also  used  when  evaluating  Equation  48.  Under 
the  assumption  that  h{x',y')  is  a  discrete  and  periodic  function,  the  Fast  Fourier 
Transform  (FFT)  can  be  applied: 

/  oo 

r«.,,„,.(Pi.P2)  =  2nGApFFT-'  ^  FFr(A"),.,„ 

V  n=l 

(49) 

where  pi  =  0,  ...,Mi  —  1,  p2  =  0,  ...,M2  —  1,  and  Mi,  M2  are  the  total  number  of 
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samples  in  the  x  and  y  directions,  respectively,  with; 


fi 


pi 


Pi 

Ax'Mi  ’ 


f‘ip2 
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,  for  Pi 
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-  1,  P2 


2  ’ 


(50) 


where  Ax'  and  Ay'  are  the  sample  intervals  in  the  x  and  y  directions,  respectively.  For 
a  derivation  of  the  spectral  component’s  conjugate  symmetry  satisfaction  requirement, 
refer  to  [59]. 

Unfortunately,  there  are  three  drawbacks  of  using  the  method  based  on  Parker’s 
work.  First,  the  assumption  of  discrete  and  finite  data  causes  biases  in  the  diagonal 
components  of  the  gradient  tensor  [61].  Second,  the  computation  point,  or  gradiome- 
ter  altitude,  must  be  held  constant.  Furthermore,  the  gradiometer  altitude  must 
also  be  above  the  highest  elevation  in  the  grid.  In  other  words,  for  AGL  type  map 
generation,  this  method  will  not  produce  reliable  results  and  a  rigorous  rectangular 
prism  method  must  instead  be  used.  While  the  drawbacks  of  this  method  have  been 
highlighted,  the  primary  advantage  is  a  rapid  reduction  in  computation  time  of  the 
gradients  maps.  Table  2  shows  the  computation  times  for  producing  full  tensor  grid- 
ded  gravity  gradient  data  via  the  rectangular  prism  and  Parker  methods.  For  more 
information  regarding  various  map  making  techniques,  refer  to  [61]  and  [62]. 


Table  2:  Gridded  Gradient  Gomputation  Time  Gomparison 


Grid 

Number  of  Points 

Rectangular  Prism 
Time 

Parker’s  Method 
Time 

3°  X  3°,  1  arc  min 

32761 

2640s 

7s 

2°  X  2°,  3  arc  sec 

5764801 

DNF* 

650s 

DNF=Did  Not  Finish  -  simulation  was  terminated  after  72  hours 
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800. 


Figure  26:  Example  Elevation  Grid. 

Using  2°  by  2°  SRTM  Digital  Terrain  Elevation  Data  (DTED)  supplied  by 
the  NGA,  discrete  elevation  grids  with  approximately  90  meter  grid  spacing  (3  arc 
second)  will  be  constructed.  Using  a  flat-earth  approximation,  the  “bottom”  of  the 
bounding  plane  is  defined  by  0  meters  MSL  which  is  assumed  to  be  located  on  the 
EGM96  geoid  surface.  Gonsequently,  the  topographical  surface  is  the  terrain  height 
of  each  grid  cell  in  meters  above  MSL.  A  small  sample  of  the  discrete  elements  used  to 
compute  gradients  due  to  terrain  is  shown  in  Figure  A.l.  For  this  method,  a  constant 
terrain  density  of  2.67g/cm^  is  assumed.  This  is  considered  by  geologists  to  be  the 
average  terrain  density.  As  previously  mentioned,  geological  surveys  have  shown  that 
actual  terrain  mass  has  varying  density.  To  mitigate  this  effect,  the  assumption  that 
the  modeled  gradient  maps  represent  truth  is  made  for  this  study. 

At  low  altitudes,  much  of  the  gradient  disturbance  signal  is  caused  by  the  terrain 
in  the  immediate  vicinity.  As  such,  relatively  small  grid  sizes  can  be  used  to  capture 
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most  of  the  terrain  effects.  Zhu  [61]  showed  that  at  an  altitude  of  10m  above  the 
maximum  altitude  of  a  relatively  rough  grid,  an  area  extent  of  a  half  degree  was 
needed  to  achieve  an  accuracy  of  lEo.  However,  as  altitude  increases,  the  area  of 
influence  on  the  signal  from  the  terrain  will  increase.  Be  that  as  it  may,  the  gradients 
from  the  terrain  are  also  falling  rapidly  as  the  altitude  increases.  For  this  study,  it  is 
assumed  that  the  grids  are  adequately  sized  for  reliable  terrain  gradient  results.  To 
minimize  gradient  errors  caused  by  grid  edge  effects,  only  the  central-most  portions 
of  the  gradient  grids  will  be  used  for  simulations. 

Gravity  Gradiometer  Modeling 

While  gravity  gradiometers  designed  for  airborne  surveys  are  inherently  com¬ 
plex,  the  three  main  drivers  of  the  signal  they  produce  are  gradient  production  rate, 
bandwidth  after  hltering,  and  noise.  Because  modeling  the  inner  workings  of  a  GGI 
are  beyond  the  scope  of  this  feasibility  study  (and  often  proprietary),  the  method 
used  to  simulate  the  signal  will  involve  manipulation  of  the  gradient  maps.  To  get 
actual  gradient  values  from  the  map,  a  table  lookup  function,  based  on  the  velocity  of 
the  GGI  host  vehicle  that  is  traveling  across  the  map,  is  performed  at  the  appropriate 
sampling  rate  using  Simulink@.  If  the  location  falls  between  grid  points,  a  spline 
interpolation  is  used  [63].  Next,  appropriate  noise  with  respect  to  the  sensor  gradient 
production  rate  will  be  added  to  the  gradient  maps  [64].  Finally,  a  low  pass  hlter 
with  similar  specihcations  to  those  used  on  gradiometer  data  will  be  implemented  on 
noisy  samples  taken  from  the  gradient  maps  to  simulate  methods  of  noise  reduction. 
In  essence,  this  hlter  will  serve  to  reject  higher  frequencies,  thus  smoothing  the  signal. 
Two  gradiometers  were  chosen  for  modeling  based  on  the  literature  review  in  Ghapter 
2. 


Noise  Generation.  As  mentioned  in  Ghapter  2,  gradiometer  manufacturers 
claim  zero-mean  gaussian  white  noise  characteristics  for  their  GGIs  over  a  certain 
bandwidth.  These  specihcations  are  given  in  terms  of  a  noise  spectral  density  (NSD). 
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The  NSD  is  the  power  of  noise  over  a  given  range  of  frequencies.  Since  the  noise 
is  white,  the  power  is  flat  over  all  applicable  frequencies.  In  order  to  determine  the 
RMS  value  of  the  white  noise  produced  by  a  GGI  before  hnal  hltering,  the  following 
relationship  is  used: 


RMS  Noise{Eo) 


2 

X 


1 

2 


Gradient  Production  Rate{Hz), 


(51) 


where  it  is  assumed  that  the  gradiometer  noise  spectral  density  is  valid  from  OHz  (or 
“DG”)  to  the  Nyquist  frequency  (defined  here  as  1/2  the  gradient  production  rate). 
Once  the  RMS  value  of  the  noise  is  found,  it  will  be  added  to  the  gradient  maps  via 
the  “normrnd”  command  in  Matlab®.  Since  the  noise  is  zero-mean,  the  RMS  values 
are  equal  to  the  standard  deviation,  a. 


Filtering.  No  matter  the  type  of  gradiometer,  all  gradient  signals  are  sent 
through  a  hnal  low  pass  hlter  (LPF)  to  reduce  uncompensated  noise  and  to  prevent 
or  reduce  aliasing  of  the  signal  [23,25].  One  drawback  to  this  type  of  hltering  is  a 
smoothing  ehect  on  the  signal.  In  other  words,  spatial  resolution  will  be  lost  for  the 
sake  of  noise  reduction.  The  most  commonly  used  LPF  in  GGI  data  noise  reduction  is 
known  as  a  Butterworth  hlter  [24,25,61]  and  its  transfer  function,  H,  is  given  by  [65]: 


H? 


(52) 


where  fc  is  the  cutoh  frequency  and  n  is  the  order  of  the  hlter. 

An  example  of  a  one  dimensional  Butterworth  hlter  with  a  normalized  cutoh 
frequency  of  0.4  and  varying  orders  is  shown  in  Figure  27.  The  frequency  spectrum 
that  is  allowed  to  pass  through  the  hlter  is  known  as  the  passband  while  the  spectrum 
that  is  cutoh  is  known  as  the  stopband.  On  an  ideal  LPF,  the  terminal  slope,  or  roll  oh. 
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CO  (normalized  to  Nyquist  frequency) 

Figure  27:  Example  Butterworth  Filter. 

between  the  passband  and  stopband  would  be  a  vertical  line  but  in  practice  a  physical 
circuit  cannot  generate  this  type  of  response.  While  Butterworth  hlters  tend  to  roll 
off  more  slowly  than  other  types  of  low  pass  hlters  (such  as  the  Chebyshev)  they  have 
very  low  ripple  characteristics.  Note  that  regardless  of  the  hlter  order,  the  magnitude 
response  of  the  Butterworth  hlter  is  always  3dB  down  at  the  cutoh  frequency.  Also 
note  the  non-constant  relationship  between  frequency  and  phase  within  the  passband 
-  one  of  the  drawbacks  of  the  Butterworth  hlter. 

For  this  research,  a  digital  7th  order  Butterworth  hlter  (Visser,  Murphy,  Lane) 
was  designed  in  Matlab(g)  and  applied  real-time  as  the  simulated  gravity  gradient 
were  traversed.  Real-time  hlter  application  marks  a  departure  from  the  considerable 
post-mission  processing  that  is  done  with  most  current  GGI  data.  For  navigation,  one 
does  not  have  the  luxury  of  extensive  post-hight  processing  -  the  hltered  signal  must 
be  available  immediately.  Real  time  application  of  the  Butterworth  hlter  as  signals  are 
being  sampled  gives  a  recursive  ehect.  That  is,  the  hlter  is  auto-regressive  in  that  it 
relies  on  previously  hltered  samples  to  compute  new  values  (i.e.  feedback).  With  some 
assumptions  regarding  the  impulse  response  of  the  hlter,  it  can  be  considered  to  give 
a  moving-average  ehect.  As  such,  this  will  inevitably  give  the  real-time  hltered  signal 
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a  lag  in  the  time  domain  as  previous  values  affect  the  most-recent  one.  Additionally, 
the  non  constant  phase  lag  of  the  Butterworth  hlter  will  cause  a  distortion  in  the 
already  lagging  hltered  signal  and  may  be  an  issue  for  a  map-matching  algorithm. 
The  hlter  function  of  the  Butterworth  LPF  has  the  following  general  form: 


N  N 

x\i)  =  akx{i  —  k)  —  bix'{i  —  1),  (53) 

k=0  1=1 


where  and  bk  are  the  hlter  coefficients,  x{i)  is  the  i-th  raw  sample,  x'{i)  is  the  i-th 
hltered  sample,  and  N  is  the  order  of  the  hlter.  It  is  evident  that  for  a  7th  order  hlter, 
it  will  require  raw  and  hltered  information  from  the  previous  7  samples.  Since  hltered 
information  is  unavailable  when  the  hlter  is  initially  applied,  it  is  assumed  that  the 
hltered  samples  are  equal  to  the  raw  samples  for  the  hrst  7  samples. 

A  key  point  when  dealing  with  hltering  on  a  moving  platform  is  that  sensed 
wavelengths  are  a  function  of  velocity^  v.  The  relationship  between  frequency  (in 
units  of  Hz),  /,  and  spatial  frequency  (in  units  of  cycles /m),  f spatial,  is  given  by: 


V 


(54) 


where  v  is  the  relative  velocity  between  the  sensor  and  the  object  being  measured 
(in  m/s).  The  corresponding  wavelength,  A  (in  units  of  meters),  is  given  by: 


A 


1 


f spatial 


(55) 


In  other  words,  the  slower  the  host  vehicle  is  moving  relative  to  the  ground,  the 
shorter  the  wavelength  (higher  frequency)  that  the  LPF  will  allow  to  be  sensed.  Given 
a  hxed  cutoff  frequency,  this  gives  the  most  resolution  and  is  why  airborne  gravity 
gradient  surveys  are  generally  flown  as  slow  as  safely  possible.  Likewise,  the  faster  the 
relative  velocity,  the  less  spatial  resolution  (again  assuming  that  the  cutoff  frequency 
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of  the  filter  is  held  constant).  To  keep  high  resolutions  at  higher  velocities,  the  obvious 
choice  is  to  increase  the  cutoff  frequency  of  the  hlter.  This  has  two  major  down  sides. 
First,  the  amount  of  noise  manifesting  itself  into  a  signal  is  generally  higher  at  high 
frequencies  and  can  quickly  render  the  signal  useless.  Second,  (assuming  that  the 
measurement  device  is  making  discrete  samples)  if  the  cutoff  frequency  is  raised  to  a 
point  where  it  becomes  close  the  the  sampling  frequency,  fs,  a  phenomenon  known 
as  aliasing  can  occur.  Aliasing  is  the  inability  to  distinguish  different  parts  of  the 
frequency  spectrum  of  a  signal  due  to  sampling  rate  restrictions.  It  is  caused  by 
having  frequency  content  of  a  signal  that  is  >  -fs-  For  example  (Figure  28),  the 
function  sin^iflu  —  0.6)t)  has  a  frequency  of  IHz.  If  it  is  sampled  at  <2Hz,  no 
unique  frequency  measurement  can  be  reconstructed  and  the  higher  frequency  signal 
can  appear  as,  or  alias  to,  a  lower  frequency.  If  the  signal  is  sampled  at  a  rate 
more  than  twice  the  maximum  frequency  content  of  the  original  signal,  the  signal  can 
be  reconstructed.  This  is  known  as  the  Nyquist  condition  and,  in  addition  to  noise 
considerations,  also  drives  GGI  filter  design.  In  practice,  it  is  desired  to  sample  the 
signal  at  approximately  5  times  the  frequency  content  of  the  signal  to  better  capture 
magnitude  and  energy  information.  It  should  be  noted  that  anti-aliasing  hlters  are 
applied  before  the  signal  is  sampled  (or  before  the  signal  is  downsampled)  and  may 
be  one  of  several  hlters  in  the  sensor.  Additionally,  the  Nyquist  frequency  is  dehned 
as  one-half  the  sampling  rate  of  a  discrete  signal  sampling  system.  In  the  case  of 
gradiometers,  raw  accelerometer  measurements  are  made  at  very  high  rates,  some 
over  lOOHz  [22].  However,  after  averaging,  internal  hltering,  scale  factor  correction 
and  other  noise  reduction  techniques,  useful  gradient  information  is  produced  on  the 
order  of  IHz  [12]. 

In  this  study,  it  is  assumed  that  the  primary  anti-aliasing  hlter  has  already 
been  applied  and  the  LPF  used  here  is  for  signal  noise  reduction.  Furthermore,  it  is 
assumed  that  host-vehicle  gradient  contributions  are  exactly  known  and  have  been 
removed  from  the  signal.  Table  3  summarizes  the  gradiometer  noise  and  bandwidth 
specihcations,  based  on  the  literature  review  in  Ghapter  2,  used  in  the  simulations. 


(a)  fs=lHz,  Nyquist  condition  not  met,  aliasing  present 


(b)  fs=2.5Hz,  Nyquist  condition  met,  no  aliasing 
Figure  28:  Illustration  of  Aliasing  and  the  Nyquist  Condition. 

GGIl  is  a  lower  noise  sensor  that  represents  the  more  optimistic  end  of  gradiometer 
performance  expected  to  be  available  within  a  decade.  Likewise,  GGI2  is  the  higher 
noise  sensor  and  represents  the  level  of  performance  that  has  already  been  demon¬ 
strated  in  tests.  Subsequently,  GGIl  and  GGI2  may  be  referred  to  as  “low  noise”  and 
“noisy”  sensors,  respectively. 


Table  3:  GGI  Specihcations 


GGI 

NSD 

fs 

RMS 

Noise 

fc 

RMS  Noise 
after  Filtering 

1 

0.223Eo/^/Hz 

IHz 

0.158AO 

0.2Hz 

O.lEo 

2 

2.23Eo/VHz 

IHz 

1.58AO 

0.2Hz 

l.OEo 

59 


Gradiometer  Model  Verification.  Since  no  real-time  filtered  gravity  gradiome- 
ter  data  is  published  in  open  literature,  a  build-up  approach  using  available  data  will 
be  used  to  validate  the  GGI  sensor  model.  First,  to  ensure  Parker’s  method  is  im¬ 
plemented  correctly,  results  from  the  method,  using  a  20th  order  expansion,  will  be 
validated  against  gradients  from  the  summation  of  the  rigorous  calculations  of  each 
rectangular  prism’s  contribution  to  the  overall  gradient.  These  rigorous  calculations 
are  a  summation  of  each  prism’s  effect  calculated  via  Equations  8-12.  A  330  x  260/cm, 
1  arc  minute  spaced  (~  1.8/cm)  grid  will  be  used  for  this  comparison.  The  fairly  large 
1  arc  minute  spacing  was  chosen  due  to  the  computational  expense  of  the  rigorous 
rectangular  prism  method.  As  a  hnal  validation  of  the  terrain  implied  gradient  map 
making  technique,  results  will  be  compared  to  those  derived  by  Zhu  and  those  cal¬ 
culated  by  Bell  Geospace  for  a  track  surveyed  by  an  Air-FTG  gradiometer  during  a 
flight  in  2004.  Zhu’s  DEM  gradient  data  is  based  on  a  numerical  integration  method 
using  a  USGS  provided  1°  x  1°,  1  arc  second  elevation  grid.  Bell  Geospace’s  ter¬ 
rain  implied  gradient  data  is  calculated  from  SRTM  data  in  the  area.  Since  actual 
data  was  unable  to  be  obtained  for  the  Zhu  and  the  Bell  Geospace  gradient  models, 
data  from  Parker’s  method  will  be  calculated  along  the  same  track,  converted  into 
the  appropriate  coordinate  frame,  plotted  and  superimposed  over  a  comparison  hgure 
originally  taken  from  Zhu  [61].  Since  EGM96  data  will  vary  depending  on  the  type  of 
interpolation  used  between  points,  results  from  the  overall  gradient  modeling  effort 
will  be  compared  to  plots  of  Bell  Geospace  Air-FTG  survey  data  in  a  manner  similar 
to  the  Parker’s  method  validation. 

For  hlter  and  noise  validation,  the  hlter’s  response  to  a  unit  impulse  input  will 
hrst  be  examined  and  compared  to  expected  values  found  using  a  method  described 
by  Rorabaugh  [66].  According  to  Rorabaugh,  the  time  at  which  the  maximum  value 
of  a  7th  order  Butterworth  hlter’s  impulse  response  occurs  should  be  approximately 
4  seconds  if  a  cutoff  frequency  of  0.2Hz  is  used.  Additionally,  the  amplitude  of  this 
peak  should  be  approximately  0.4.  To  validate  the  noise  in  the  signal,  the  mean  and 
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standard  deviation  of  the  “normrnd”  generated  noise  will  be  examined  before  and 
after  filtering  is  applied  and  shonld  match  valnes  listed  in  Table  3. 

Navigation 

To  determine  the  feasibility  of  nsing  a  GGI-based  map-matching  navigation 
system,  some  assnmptions  abont  the  system  will  be  made  in  order  to  better  classify 
whether  the  modeled  GGI  signal  exceeds  threshold  reqnirements.  It  is  assnmed  that 
the  GGI  will  be  nsed  in  a  navigation  scheme  similar  to  those  presented  by  Jircitano  [44] 
and  Richeson  [51].  That  is,  an  INS  will  be  nsed  as  the  primary  navigator  and  it’s 
error  will  be  bonnd  by  npdates  from  a  GGI  signal  to  map  matching  algorithm.  Fignre 
29  shows  a  schematic  of  a  generic  GGI/map  matching  npdated  INS. 


Fignre  29:  GGI-Aided  Passive  Navigation  System  Flowchart. 

Test  Conditions.  Two  track  areas,  one  with  with  rapid  terrain  changes  and 
and  one  with  relatively  low  levels  of  terrain  changes  (i.e.  monntainons  versns  flat 
terrain)  were  selected  for  the  navigation  feasibility  simnlations.  Track  1  (Fignre  30a) 
is  a  west-east  track  beginning  jnst  sontheast  of  Monterey,  GA,  (36.25°A^,  121.5°IF) 
in  a  region  of  relatively  monntainons  terrain.  This  area  was  selected  in  part  becanse 
Rice  University  contracted  Bell  Geospace  to  fly  a  gravity  gradient  snrvey  in  that  area 
in  2004  [67].  As  snch,  airborne  gradiometer  data  is  available  for  that  area  and  may 
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Longitude 


(a)  Track  1,  rough  terrain. 


(b)  Track  2,  smooth  terrain. 


Figure  30:  Test  Tracks  (note  contour  scale  differences). 
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be  valuable  for  future  navigation  research  efforts  given  its  vicinity  to  the  Air  Force 
Flight  Test  Center  at  Edwards  Air  Force  Base,  CA.  Track  2  (Figure  30b)  is  also  a 
west-east  track  located  in  western  Tennessee  beginning  slightly  north  of  Memphis. 
This  track,  starting  at  36°N,  89.75°fF,  was  selected  as  it  represents  one  of  the  largest 
areas  of  relatively  flat  terrain  in  the  United  States.  Note  the  different  scales  in  Figure 
30.  The  along-track  terrain  statistics  for  ~135/cm  tracks  in  each  area  are  shown  in 
Table  4.  Should  an  area  of  smooth  terrain  be  desired  for  flight  test,  the  El  Centro 
complex  near  the  Salton  sea  could  be  useful  as  it  contains  relatively  large  stretches  of 
flat  terrain. 


Table  4:  Along-Track  Terrain  Statistics 


Track 

Min 

Max 

Mean 

Std.  Dev. 

1,  Rough 

78m 

1118m 

403m 

254m 

2,  Smooth 

72m 

146m 

105m 

24m 

For  each  track,  straight  and  level  runs  with  altitudes  ranging  from  1000-20000m 
height  above  average  terrain  (HAAT)  and  velocities  from  50-1200  m/s  will  be  simu¬ 
lated  for  100  seconds  beginning  once  the  hlter  has  had  time  to  produce  meaningful 
results.  These  test  points  were  chosen  in  order  to  map  a  wider  range  of  the  flight 
envelope  than  has  been  done  in  previous  studies  [44,51].  Note  that  most  conventional 
military  aircraft  generally  operate  in  the  5-lOkm  altitude  region  at  speeds  of  100- 
300m/s.  However,  some  surveillance  aircraft  operate  at  higher  altitudes  and  loiter  at 
slower  velocities,  hence  the  inclusion  of  20km  altitudes  and  velocities  down  to  50m/s. 
Additionally,  the  GGI  signal  sensed  on  board  a  theoretical  high  speed  vehicle  may 
be  useful  for  future  studies,  thus  600-1200m/s  velocities  are  also  included.  A  time  of 
100s  was  chosen  under  the  assumption  that  the  map-matching  algorithm  will  provide 
1  second  updates  to  the  INS.  Thus,  general  trends  in  the  signal’s  usefulness  should  be 
able  to  be  seen  over  the  course  of  100s.  Table  5  summarizes  the  planned  test  runs. 
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Table  5:  Test  Matrix. 


Runs 

Altitude 
(m  HAAT) 

Velocity  (m/s) 

Remarks 

1  -  6 

1000 

50,100,150,300,600,1200 

GGI  1,  Rough  Terrain 

7-12 

1000 

50,100,150,300,600,1200 

GGI  1,  Smooth  Terrain 

13  -  18 

1000 

50,100,150,300,600,1200 

GGI  2,  Rough  Terrain 

19  -  24 

1000 

50,100,150,300,600,1200 

GGI  2,  Smooth  Terrain 

25  -  30 

2500 

50,100,150,300,600,1200 

GGI  1,  Rough  Terrain 

31  -  36 

2500 

50,100,150,300,600,1200 

GGI  1,  Smooth  Terrain 

37-42 

2500 

50,100,150,300,600,1200 

GGI  2,  Rough  Terrain 

43  -  48 

2500 

50,100,150,300,600,1200 

GGI  2,  Smooth  Terrain 

49  -  54 

5000 

50,100,150,300,600,1200 

GGI  1,  Rough  Terrain 

55  -  60 

5000 

50,100,150,300,600,1200 

GGI  1,  Smooth  Terrain 

61  -  66 

5000 

50,100,150,300,600,1200 

GGI  2,  Rough  Terrain 

67-72 

5000 

50,100,150,300,600,1200 

GGI  2,  Smooth  Terrain 

73  -  78 

10000 

50,100,150,300,600,1200 

GGI  1,  Rough  Terrain 

79  -  84 

10000 

50,100,150,300,600,1200 

GGI  1,  Smooth  Terrain 

85  -  90 

10000 

50,100,150,300,600,1200 

GGI  2,  Rough  Terrain 

91  -  96 

10000 

50,100,150,300,600,1200 

GGI  2,  Smooth  Terrain 

97  -  102 

20000 

50,100,150,300,600,1200 

GGI  1,  Rough  Terrain 

103  -  108 

20000 

50,100,150,300,600,1200 

GGI  1,  Smooth  Terrain 

109  -  114 

20000 

50,100,150,300,600,1200 

GGI  2,  Rough  Terrain 

115  -  120 

20000 

50,100,150,300,600,1200 

GGI  2,  Smooth  Terrain 

121 

5000 

150 

GGI  1,  Rough,  Form  w/KG-10 

122 

5000 

150 

GGI  1,  Smooth,  Form  w/KG-10 

123 

5000 

150 

GGI  2,  Rough,  Form  w/KG-10 

124 

5000 

150 

GGI  2,  Smooth,  Form  w/KG-10 

125 

10000 

150 

GGI  1,  Rough,  Form  w/KG-10 

126 

10000 

150 

GGI  1,  Smooth,  Form  w/KG-10 

127 

10000 

150 

GGI  2,  Rough,  Form  w/KG-10 

128 

10000 

150 

GGI  2,  Smooth,  Form  w/KG-10 

The  Tanker  Effect.  Since  aerial  refueling  (Figure  31)  is  an  integral  part  of 
military  force  projection,  it  is  desired  to  know  if  and  by  what  amount  the  presence 
of  a  large  tanker  aircraft  in  the  vicinity  of  the  GGI  carrying  aircraft  will  corrupt  the 
GGI  signal.  While  it  typically  takes  approximately  hve  minutes  to  refuel  a  fighter 
size  aircraft,  these  aircraft  will  often  stay  in  formation  with  the  tanker  for  extended 
periods  of  time,  especially  during  transit  to  forward  locations. 
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Figure  31:  A  KC-10  Offloading  Fuel  to  an  F-22  Raptor. 


To  model  the  gravitational  gradient  effects  of  a  tanker  aircraft,  a  simple  rect¬ 
angular  prism  model  of  a  McDonnell  Douglas  KC-10  is  created  based  on  dimensions 
derived  from  Figure  32,  taken  from  [68].  It  should  be  noted  that  the  volume  calcula¬ 
tion  is  only  an  approximation  based  on  the  relative  dimensions  of  the  KC-10.  Once 
the  approximate  volume  is  known,  a  typical  heavy  weight  for  a  KC-10  in  flight  is 
used  to  calculate  the  average  density.  Finally,  to  obtain  a  representative  rectangular 
prism,  the  dimensions  required  to  obtain  the  approximate  volume  were  best  fitted  to 
the  KC-10  as  shown  in  Figure  33.  All  estimated  properties  are  listed  in  Table  6.  While 
a  higher  resolution  model  of  the  tanker  could  have  been  created  via  more  rectangular 
prisms,  the  primary  objective  here  is  only  to  determine  if  a  large  object,  located  next 
to  and  flying  at  the  same  velocity  as  the  sensor,  having  a  density  distribution  and 
size  that  represents  a  large  aircraft,  will  affect  the  signal.  For  this  study,  the  GGI  will 
nominally  be  located  30m  laterally  and  co-altitude  with  respect  to  the  center  of  the 
rectangular  prism,  but  allowed  to  vary  approximately  ±5m  in  all  directions  to  sim- 
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Table  6:  KC-10  Parameters 


Parameter 

Value 

Total  Mass 

240,000A:c/ 

Approx.  Volume 

1816m^ 

Approx.  Density 

132^ 

mr 

Length  of  Equiv.  Rect.  Prism 

48.8m 

Width  of  Equiv.  Rect.  Prism 

6.1m 

Height  of  Equiv.  Rect.  Prism 

6.1m 

ulate  relative  motion  between  the  two  aircraft  typically  seen  during  formation  flight. 
Gradients  will  be  calculated  at  the  8  points  which  dehne  the  boundary  of  the  GGI 
aircraft’s  deviation  from  the  nominal  position.  The  mean  and  standard  deviation  of 
these  values  will  be  computed  and  the  “normrnd”  command  will  be  used  to  convert 
these  gradients  into  effective  white  noise  which  will  then  be  added  to  the  gradient 
maps.  Next,  the  hlter  will  be  applied  as  the  maps  are  traversed.  The  hltered  signals 
with  and  without  the  tanker  effects  will  then  be  plotted  for  comparison.  Addition¬ 
ally,  these  runs  will  only  be  accomplished  at  typical  refueling/cruise  altitudes  and 
velocities  (Table  5). 


Figure  32:  KG-10  Dimensions. 
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6.1ni 


P  =  132  kg/m^ 


6.1m 


Figure  33:  KC-10  Transformed  into  a  Rectangular  Prism. 
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Gradiometer  Signal  Metrics  for  Useful  Navigation.  For  a  signal  to  be  usable 
in  a  map  matching  utility  sense,  it  must  vary  in  time  at  an  adequate  rate  while 
maintaining  good  signal  to  noise  characteristics.  If  the  signal  doesn’t  vary  much  with 
time,  or  there  is  excessive  noise,  matching  it  to  a  map  would  prove  to  be  an  exercise  in 
futility  and  tell  very  little  about  relative  position  on  the  map.  In  order  to  determine  if 
the  simulated  GGI  signal  is  useful  for  navigation,  several  signal  examination  methods 
will  be  employed.  First,  the  signal  will  be  analyzed  graphically  to  determine  if  the 
along-track  signal  changes  are  outside  of  the  noise  level  from  the  gradiometer.  To 
do  this,  the  standard  deviation  of  the  noise  will  be  plotted  in  relation  to  the  mean 
value  of  the  true  gradient  over  the  run  time.  Additionally,  the  overall  uniqueness 
and  clarity  of  the  gradient  contour  outlined  by  the  GGI  will  be  examined.  If  less 
than  approximately  5%  of  the  signal  is  within  the  standard  deviation  of  the  noise 
and  the  contour  is  clearly  defined,  the  signal  is  considered  “excellent”.  That  is,  the 
signal  defines  the  contour  well  and  the  contour  varies  in  time  to  provide  uniqueness. 
Should  approximately  5-20%  of  the  signal  fall  outside  the  noise  standard  deviation 
and  the  contour  remain  easily  discernable,  the  signal  will  be  deemed  “useful”.  If 
approximately  20-50%  of  the  signal  is  outside  the  standard  deviation  and  contour 
uniqueness  begins  to  drop  but  is  still  somewhat  discernable,  the  signal  is  considered 
“marginally  useful”.  Finally,  if  greater  than  50%  of  the  signal  is  within  the  noise 
standard  deviation,  the  signal  is  considered  “unusable”.  That  is,  the  signal  is  too 
noisy  to  uniquely  define  the  contour.  These  metrics  bear  no  qualitative  data  as  to 
expected  navigation  performance  but  are  designed  to  give  an  idea  of  the  relative  map¬ 
matching  usefulness  of  the  signal  at  different  flight  conditions.  While  a  signal  to  noise 
ratio  (SNR)  examination  was  considered,  a  strong  SNR  does  not  guarantee  that  the 
signal  varies  in  time. 

Next,  a  comparison  of  signal  levels  obtained  from  GGIl  will  be  measured  against 
those  previously  proven  useful  for  navigation  by  Richeson.  GGIl  is  chosen  because 
Richeson  found  that  a  gradiometer  with  lEo  of  RMS  error  (i.e.  GGI2)  provided 
no  gains  over  unaided  navigation  grade  IMUs  when  used  with  his  map-matching 


Table  7:  Richeson’s  “GGI  Survey”  simulation  parameters 


Parameter 

Value 

Track 

45.0° A,  113.0°  -  100.3°IU 

Velocity 

40m/s 

Altitude 

100m 

Gradient  Model 

EGM96  based 

dTij/dt{RMS) 

0.05Eo/s 

algorithm.  However,  a  GGI  with  0.1i7o  of  RMS  noise  did  provide  some  benefit  to  the 
INS,  though  not  to  GPS  levels.  Though  terrain  effects  were  neglected  and  only  EGM96 
derived  gradient  maps  were  produced,  Richeson’s  research  indirectly  classified  the 
signal  threshold  required  for  an  improvement  in  navigation  over  an  unaided  navigation 
grade  IMU,  assuming  a  gradiometer  with  O.lEo  of  noise  [51].  The  metrics  which  define 
this  threshold  are  obtained  by  examination  of  the  average  along  track  signal  rate  of 
change  for  each  gradient,  before  noise  is  added. 

Figure  34  shows  the  differences  in  navigation  results  from  Richeson’s  research 
using  a  O.lEo  GGI  and  a  GPS  in  combination  with  navigation  grade  IMUs.  Given 
enough  time,  the  GGI  aided  navigation  system  bounds  the  INS  error  to  10-20m  in 
the  north,  east,  and  down  directions  whereas  the  GPS  bounds  the  error  of  the  system 
to  approximately  0.01m.  In  other  words,  the  GPS  aided  system  provides  results 
that  are  roughly  3  orders  of  magnitude  better  that  the  O.lEo  GGI  aided  system. 
Since  a  gradiometer  specific  map-matching  algorithm  is  still  under  development  at 
AFIT,  the  assumption  is  made  that  order  of  magnitude  differences  in  the  noise  free 
signal  time  rate  of  change  will  correspond  to  order  of  magnitude  changes  in  navigation 
performance  (i.e.  the  more  unique  the  GGI  signal  is,  the  more  likely  it  will  be  correctly 
correlated  to  a  map  and  the  better  the  navigation  performance  will  be).  Based  on 
Richeson’s  results,  a  signal  having  a  0.05Eo/s  RMS  rate  of  change  and  a  GGI  with 
O.lEo  of  noise  will  bound  host- vehicle  position  error  to  roughly  10-20m.  Using  this 
metric.  Table  8  was  constructed. 

It  must  be  stressed  that  this  only  holds  under  the  assumption  that  the  gradient 
maps  are  truth.  If  the  maps  inherently  have  error  in  them  (i.e.  if  GPS  was  used  for 
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(a)  Navigation  results  for  a  O.lEo  GGI  map-matching  system  pro¬ 
viding  IHz  updates  to  the  INS. 
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(b)  Navigation  results  for  a  GPS  providing  IHz  updates  to  the 
INS. 


Figure  34:  Richeson’s  Navigation  Results:  O.lEo  GGI  Map-Matching  System  vs  GPS, 
taken  from  [51]. 
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Table  8:  O.lEo  GGI  Signal  Glassification  Metrics 


Navigation  Usefulness 
(Anticipated  Error  Bounds) 

Signal  Time  Rate  of  Change 
RMS,  noise  free 

Unusable 

<  0.05Eo/s 

Marginally  Useful  (lO^m  MRSE) 

0.05  -  0.5Eo/s 

Useful  MRSE) 

0.5-5Eo/s 

Excellent  (10“lm  MRSE) 

>  5Eo/s 

positioning  during  a  survey),  navigation  performance,  in  a  best  case  sense,  will  be 
limited  to  the  amount  of  error  present  in  the  maps. 


Terrain  Avoidance 

The  use  of  a  gravity  gradiometer  as  a  terrain  avoidance  warning  enhancement 
is  a  distinct  challenge.  According  to  Gleason  [23],  the  three  parameters  which  dictate 
feasibility  of  the  GGI  in  the  role  of  terrain  avoidance  are  the  cutoff  frequency  of  the 
hlter  used  to  suppress  uncompensated  error  sources,  the  gradient  production  rate, 
and  the  hnal  gradiometer  noise  level.  This  is  compounded  by  the  fact  that  many 
aircraft  which  maneuver  at  very  low  altitudes  and  could  beneht  from  a  passive  terrain 
avoidance  system  typically  fly  at  velocities  often  on  the  order  of  several  hundred  meters 
per  second.  Assuming  the  GGI  bandwidth  is  fixed,  the  increase  in  speed  serves  to 
limit  the  shortest  sensed  wavelengths.  Table  9  shows  velocity  versus  minimum  sensed 
wavelength  with  a  gradiometer  cutoff  frequency  of  0.2Hz  using  the  relationship  shown 
in  Equation  55.  Since  the  terrain  makes  up  the  higher  frequency  end  of  the  overall 


Table  9:  Velocity  vs.  Minimum  Wavelength,  fc=0.2Hz 


Velocity  (m/s) 

Minimum  Sensed  Wavelength  (m) 

50 

250 

100 

500 

150 

750 

300 

1500 

600 

3000 

1200 

6000 

gradient  spectrum,  hltering  high  frequencies  out  in  the  name  of  noise  reduction  may 
prove  devastating  for  terrain  avoidance.  Gompounding  the  loss  of  spatial  frequency 
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information  is  the  fact  that  locating  the  mass  anomaly  corresponding  to  a  gradient 
is  an  inverse  problem.  The  following  question  must  be  posed:  If  a  gradient  change 
is  sensed,  is  the  change  from  a  large  mountain  in  the  distance  or  from  a  small  tower 
that  the  aircraft  in  question  is  about  to  impact? 

Method  1.  For  this  study,  the  assumption  is  made  that  no  prior  positional 
information  is  known  and  a  GGI  with  a  IHz  gradient  production  rate  and  0.2Hz 
cutoff  is  being  used  as  the  sole  device  in  an  attempt  to  provide  a  consistent  terrain 
avoidance  warning.  In  other  words,  besides  GGI-provided  information,  the  user  has 
little  situational  awareness.  Five  runs,  using  a  noise  free  version  of  the  aforementioned 
gradiometer,  will  be  flown  from  west  to  east  over  a  perfectly  flat  surface  on  which 
obstacles  of  varying  size  will  be  placed.  The  obstacles  will  be  cubic  with  dimensions 
of  25,  50,  100,  250,  and  500m,  all  having  a  density  of  2.Q7g/cm?.  The  aircraft  will  be 
flown  straight  and  level  at  50m/s  on  a  plane  10m  below  the  top  of  each  object.  The 
no  noise,  relatively  slow  velocity  characteristics  were  chosen  to  make  this  a  best-case 
scenario.  If  feasibility  is  not  demonstrated  for  this  case,  then  the  addition  of  noise 
and  higher  velocities  will  only  exacerbate  the  situation.  Figure  35  shows  the  overall 
scenario  setup.  The  geometry  will  make  the  T^x  gradient  of  utmost  interest  since  it 
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defines  the  east-west  edges  of  the  anomaly.  To  potentially  connter  the  inverse  natnre 

of  determining  the  terrain  distance  from  the  gradient,  the  time  rate  of  change  of  the 

dT 

Txx  gradient  will  also  be  examined.  A  threshold  value  of  T^z  or  — ^  corresponding  to 

at 

imminent  obstacle  impact  will  be  sought.  For  this  study,  imminent  impact  is  dehned 
as  1.5s  to  impact  [13]. 

Method  2.  The  premise  for  this  potion  of  the  study  is  that  the  navigation 
system  provides  adequate  latitude,  longitude,  and  altitude  information  and  that  this 
position  information  is  used  to  perform  lookups  on  terrain  elevation  databases  and 
gravity  gradient  maps  stored  onboard  the  aircraft.  Thus,  it  is  assumed  that  if  a 
correct  terrain  elevation  database  is  used,  impact  with  the  modeled  terrain  can  be 
prevented  (or  at  least  predicted).  The  GGI’s  role  then  becomes  to  warn  the  navigation 
system,  and  ultimately  the  operator,  of  unmodeled  terrain  anomalies  which  may  be 
encountered  along  the  flight  path.  Some  examples  include  communications  towers, 
water  tanks,  and  other  “pop-up”  structures  that  are  often  constructed  before  terrain 
databases  can  be  updated.  Though  these  structures  should  be  listed  in  the  NOTAMs, 
such  information  isn’t  always  available  in  hostile  areas.  To  investigate  the  feasibility  of 
the  GGI’s  ability  to  properly  predict  an  potential  impact  with  an  un-modeled  object, 
a  simple  case  study  using  the  world’s  tallest  water  tower,  listed  at  218  feet  tall  with 
approximately  500,000  gallons  of  water,  will  be  executed  (see  Figure  36)  [69].  Using 


Figure  36:  World’s  Tallest  Water  Tower 
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an  assumed  weight  of  1  US  Gallon  of  water  of  3.785kg  and  1  USG  water  =  3785. 4cm^, 

the  density  of  water  was  calculated  to  be  approximately  l.Qg/cm?.  An  average  density 

of  3.9g/cm^  (half  that  of  mild  steel)  was  used  for  the  support  prism.  The  water  tower 

was  chosen  because  it  represents  the  larger  end  of  the  spectrum  of  un-modeled  pop-up 

structures.  If  the  gradiometer  can  predict  an  impact  with  it,  smaller  objects  will  be 

tested.  The  tower  will  be  modeled  with  two  rectangular  prisms,  one  to  simulate  the 

support  and  one  to  simulate  the  tank.  Before  the  simulation,  a  frequency  domain 

analysis  will  be  conducted  on  the  anomaly  to  determine  its  signal  structure.  Then, 

runs  will  be  flown  from  west  to  east  over  a  perfectly  flat  surface  (excluding  the  tower). 

This  geometry  will  again  make  the  T^x  gradient  of  utmost  interest.  As  before,  the 

time  rate  of  change  of  the  Txx  gradient  will  be  examined.  A  threshold  value  of  T^z 
dT 

or  — —  corresponding  to  imminent  object  impact  (1.5s)  will  be  sought.  Figure  37 
at 

shows  the  overall  test  setup. 


12.4m 


Figure  37:  Modeling  the  Unmodeled 
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IV.  Results  and  Analysis 


Overview 

This  chapter  presents  the  overall  results  and  analysis  obtained  from  the  methods 
described  in  Chapter  III.  The  overall  goal  is  to  present  data  that  suggests  or  disproves 
the  feasibility  of  using  a  gravity  gradiometer  in  the  roles  of  passive  aircraft  naviga¬ 
tion  and  terrain  avoidance.  First,  multi-step  validation  results  will  be  presented  to 
show  that  the  gradiometer  models  used  in  this  study  are  producing  signals  that  rep¬ 
resent  real-time  measurements  of  gravity  gradients  produced  by  the  earth.  Next,  the 
signals’  navigation  usefulness  over  different  terrain  variety,  altitude  and  velocity  will 
be  evaluated  and  summarized  using  metrics  developed  for  this  study.  Additionally, 
the  signal  effects  of  flight  in  the  vicinity  of  a  large  tanker  aircraft  will  be  examined. 
Finally,  the  results  of  several  GGI-based  terrain  avoidance  scenarios  will  be  presented 
and  discussed. 

Model  Validation 

The  model  validation  effort  begins  with  an  examination  of  results  from  the 
gravity  gradient  map  making  process.  Figure  38  illustrates  the  differences  between 
the  gradients  from  the  rigorous  rectangular  prism  summation  and  from  Parker’s 
method.  Clearly  shown  is  the  bias  in  the  diagonal  components  of  the  gradient  tensor 
(Txx,  Tyy,  Tzz).  While  these  biases  are  relatively  large,  they  do  not  affect  the  overall 
shape,  or  uniqueness,  of  the  diagonal  gradients  as  a  function  of  distance  traveled.  As 
such,  these  biases  are  deemed  acceptable  within  the  scope  of  this  research.  Also  shown 
is  the  excellent  correlation  of  off  diagonal  terms  generated  via  the  two  methods.  It 
should  be  noted  that  the  hrst  and  last  hfth  of  the  original  gradient  grids  were  excluded 
due  to  edge  effects  that  manifested  during  map  generation.  All  gradient  maps  used 
in  the  simulations  were  corrected  for  edge  effects. 
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Figure  38:  Model  Validation  -  Parker’s  Method  vs  Rigorous  Rectangular  Prisms. 


With  a  fast  method  of  gradient  calculation  successfully  implemented,  a  compari¬ 
son  of  results  from  this  technique  are  compared  to  those  used  by  geophysicists.  Figure 
39  shows  the  final  terrain  generated  gradient  map  verification.  The  data  illustrates 
terrain  implied  gradients  that  were  calculated  along  a  track  surveyed  by  an  airborne 
gradiometer  during  a  2004  Bell  Geospace  flight.  Note  that  Fj^  is  interchangeable  with 
Tij.  It  must  be  stressed  that  since  data  was  superimposed  onto  a  pre-existing  plot, 
this  is  a  qualitative  comparison  to  ascertain  if  Parker’s  method  is  working  properly. 
Based  on  the  plots,  is  clear  that  Parker’s  method  has  been  implemented  successfully 
and  that  the  terrain  implied  gradient  maps  generated  in  this  study  compare  to  those 
generated  by  geophysicists.  Also  evident  are  the  biases  present  on  the  diagonal  com¬ 
ponents  of  the  gradient  tensor.  The  slight  mismatch  in  the  off-diagonal  components 
is  likely  due  to  the  fact  that  the  exact  track  coordinates  were  not  known  and  were 
approximated  by  visual  examination  of  a  map. 
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Figure  39:  Terrain  Effect  Modeling  verification:  Bell  Geospace  vs.  Zliu’s  Numerical 
Integration  vs.  Parker’s  Method,  taken  from  [61]  and  modified. 


Figure  40  compares  the  complete  gravity  gradient  model  to  actual  values  ob¬ 
tained  by  a  Bell  Geospace  Air-FTG  gradiometer  [61,67].  While  the  gradients  are  not 
an  exact  match,  the  trends  are  clearly  predicted  by  the  model.  The  larger  differences 
in  this  hgure  are  likely  caused  by  the  gradiometer  sensing  density  anomalies  within 
the  terrain  and  geology  which  are  not  accounted  for  by  the  model.  As  before,  the 
exact  track  coordinates  were  unknown  -  also  likely  contributing  to  the  differences.  It 
must  be  stressed  that  this  hgure  illustrates  units  of  10“®  and  any  incorrect  coding  or 
other  error  is  likely  to  skew  modeled  values  considerably  more  than  shown.  Based 
on  these  results,  it  is  concluded  that  the  model  is  producing  gradients  that  represent 
realistic  values  produced  by  the  earth. 
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Figure  40:  Model  Validation  -  Modeled  Gradients  vs.  Air-FTG  Data,  taken  from  [61] 
and  modified. 

Figure  41  shows  the  frequency  response  of  the  7th  order  Butterworth  hlter  used 
in  this  study.  The  gradiometers  to  which  this  hlter  was  applied  have  a  gradient 
production  frequency  of  l.OHz,  thus  the  Nyquist  frequency  is  one-half  that  (0.5Hz). 
As  stated  in  the  gradiometer  specihcations  (Table  3),  the  LPF  cutoff  frequency  is 
0.2Hz  and,  when  normalized  to  the  Nyquist  frequency,  becomes  0.4,  as  shown  in  the 
hgure.  Figure  42  shows  the  impulse  response  of  the  hlter.  As  expected,  the  peak 
value  is  at  4  samples,  which,  with  a  sampling  rate  of  IHz,  corresponds  to  4s.  Also 
noted  is  the  fact  that  the  amplitude  matches  the  predicted  value  of  0.4  as  well.  Based 
on  this  simple  study,  it  is  concluded  that  the  hlter  used  in  this  model  was  properly 
implemented. 
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Figure  41:  Filter  validation  -  Frequency  response  of  a  7tli  order  Butterwortli  filter 
with  fc=0.2Hz  and  fs=lHz. 


Figure  42:  Filter  Validation  -  Impulse  Response  of  a  7th  Order  Digital  Butterworth 
Filter  with  fc=0.2Hz  and  fs=lHz. 

Figure  43  shows  the  mean  and  standard  deviation  of  100  samples  of  100  ran¬ 
dom  numbers  generated  using  the  “normrnd”  command  in  an  attempt  to  validate 
the  simulated  white  noise  before  and  after  the  hltering  used  for  each  GGI.  Glearly, 
the  mean  noise  of  both  GGIs’  pre  and  post-filtered  signals  are  approximately  zero,  as 
anticipated.  Recall  that  GGIl  represents  a  relatively  low  noise  sensor  projected  to  be 
available  in  10  years  and  that  GGI2  represent  a  more  noisy  sensor  that  is  currently 
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in  flight  test.  The  standard  deviation  of  GGIl’s  pre-hltered  noise  is  approximately 
0.158i?o  while  the  output  noise  is  approximately  Q.lEo.  Similarly,  the  standard  de¬ 
viation  of  GGI2’s  pre-hltered  noise  is  approximately  1.58Eo  while  the  output  noise  is 
approximately  lEo.  All  values  are  expected  and  indicate  that  the  noise  generation 
process  and  hlter  are  working  correctly. 


(a)  GGI  1  noise. 
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(b)  GGI  2  noise. 

Figure  43:  GGI  Noise  Validation. 
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Aircraft  Navigation  Feasibility  Results 


In  an  effort  to  condense  the  results,  only  the  “zz”  (or  down-down)  component 
of  the  gravitational  disturbance  gradients  will  be  examined  in  this  section.  This 
component  was  chosen  because  it  is  the  most  intuitive  as  it  generally  relates  to  the 
shape  of  the  terrain.  The  remaining  independent  components  {T^x,  Txy,  Tyy  and  Tyf) 
behave  similarly  and  will  offer  additional  uniqueness  to  the  overall  map-matching 
effort.  Though  they  are  omitted  here,  the  data  is  available  for  future  research.  Plots 
of  Tzz  for  all  navigation  runs  (1-128)  are  located  in  Appendix  C. 


Navigation  Feasibility  via  Qualitative  Signal  Analysis.  Tables  10  and  11 
summarize  the  signal  characteristics  derived  from  simulation  results  using  GGIl,  the 
lower  noise  gradiometer. 


Table  10:  GGIl  Signal  Glassihcation  Results  -  Rough  Terrain 


50m/s 

lOOm/s 

150m/s 

300m/s 

600m/s 

1200m/ s 

1000m  HA  AT 

E 

E 

E 

E 

U 

X 

2500m  HA  AT 

E 

E 

E 

E 

E 

U 

5000m  HAAT 

E 

E 

E 

E 

E 

E 

10000m  HAAT 

E 

E 

E 

E 

E 

E 

20000m  HAAT 

U 

E 

E 

E 

E 

E 

HAAT=Height  Above  Average  Terrain 


Table  11:  GGIl  Signal  Glassification  Results  -  Smooth  Terrain 


50m/s 

lOOm/s 

150m/s 

300m/s 

600m/s 

1200m/ s 

1000m  HAAT 

E 

E 

E 

E 

U 

MU 

2500m  HAAT 

MU 

MU 

U 

E 

E 

E 

5000m  HAAT 

X 

X 

X 

X 

U 

U 

10000m  HAAT 

X 

X 

X 

X 

u 

u 

20000m  HAAT 

X 

X 

X 

X 

MU 

u 

E=Exce’ 

[lent,  U=Useful,  MU 

= Marginal 

[y  Useful,  X=Unusable 

Figure  44  graphically  represents  the  data  from  Tables  10  and  11.  Glearly  shown 
is  the  excellent  signal  produced  by  GGIl  at  nearly  all  tested  flight  conditions  when  it 
is  flown  over  rough  terrain  (track  1).  Note  that  the  low  altitude,  high  velocity  signal 
degradation  is  caused  by  the  LPF  and  will  be  later  discussed.  The  signal  quickly 
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(b)  GGIl  signal  classification  results,  smooth  terrain. 


Figure  44:  Low  Noise  GGI  Signal  Glassification  Summary  (GGIl). 
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degrades  when  the  instrument  is  flown  over  smooth  terrain  (track  2),  as  shown  in 
Figure  44b.  Note  the  detrimental  effects  of  slower  velocities  and  increasing  altitude. 
To  put  these  results  into  real-world  terms,  most  conventional  military  aircraft  would 
fall  into  the  marginally  useful  to  unusable  regimes  if  using  this  GGI  over  smooth 
terrain.  Though  it  can  be  argued  that  un-modeled  density  variations  are  certain  to 
exist  within  and  below  the  terrain  and  would  help  to  give  uniqueness  to  the  signal, 
this  case  represents  a  worst  case  scenario  where  the  aircraft  may  be  flying  over  a  deep 
ocean  or  other  areas  of  sparse  terrain  and  relatively  constant  subterranean  density. 
Glearly,  a  more  sensitive  gradiometer  will  be  required  if  accurate  navigation  is  to  be 
maintained  in  these  conditions. 

Tables  12  and  13  summarize  the  signal  characteristics  derived  from  simulation 
results  using  GGI2,  the  noisier  gradiometer. 


Table  12:  GGI2  Signal  Glassification  Results  -  Rough  Terrain 


50m/s 

lOOm/s 

150m/s 

300m/s 

600m/s 

1200m/ s 

1000m  HA  AT 

E 

E 

E 

E 

U 

X 

2500m  HA  AT 

E 

E 

E 

E 

E 

U 

5000m  HAAT 

U 

U 

U 

U 

E 

E 

10000m  HAAT 

X 

MU 

U 

U 

E 

E 

20000m  HAAT 

X 

X 

X 

MU 

MU 

U 

HAAT=Height  Above  Average  Terrain 


Table  13:  GGI2  Signal  Glassification  Results  -  Smooth  Terrain 


50m/s 

lOOm/s 

150m/s 

300m/s 

600m/s 

1200m/ s 

1000m  HAAT 

X 

X 

X 

X 

X 

MU 

2500m  HAAT 

X 

X 

X 

X 

X 

X 

5000m  HAAT 

X 

X 

X 

X 

X 

X 

10000m  HAAT 

X 

X 

X 

X 

X 

X 

20000m  HAAT 

X 

X 

X 

X 

X 

X 

E=Exce’ 

[lent,  U=Useful,  MU 

= Marginal 

[y  Useful,  X=Unusable 

Figure  45  graphically  represents  the  data  from  Tables  12  and  13.  Glearly  shown 
is  the  overall  degradation  of  the  signal  produced  by  GGI2.  Over  rough  terrain,  the 
signal  usefulness  has  some  similarities  to  GGIl’s  performance  over  smooth  terrain. 
As  expected,  the  increase  in  GGI  noise  has  the  same  effect  as  increasing  altitude  or 
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(a)  GGI2  signal  classification  results,  rough  terrain. 
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(b)  GGI2  signal  classification  results,  smooth  terrain. 


Figure  45:  Noisier  GGI  Signal  Glassification  Summary  (GGI2). 
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decreasing  velocity.  In  real-world  terms,  most  conventional  military  aircraft  would 
fall  into  the  useful  to  marginally  useful  regimes  if  using  this  GGI  over  rough  terrain. 
When  GGI2  is  used  over  smooth  terrain,  the  noise  level  causes  a  total  loss  in  signal 
usefulness  for  all  flight  conditions  except  for  a  theoretical  low  altitude,  hypersonic 
case.  These  results  show  that  a  GGI  with  l.QEo  of  noise  will  not  be  adequate  for 
all  flight  and  terrain  conditions,  in  agreement  with  previous  works  [51].  For  added 
physical  insight,  examples  of  the  LPF  effect  and  signals  at  each  level  of  classihcation 
are  now  presented. 


Figure  46:  Low  Pass  Filter  Effect,  1200m/s,  1000m,  Rough  Terrain. 

Figure  46  shows  the  unique  case  where  high  velocity  combined  with  rapid 
changes  in  gradient  (high  frequency)  actually  cause  a  loss  in  signal  usefulness.  This 
run  was  completed  at  1000m  HAAT  and  1200m/s  (run  6).  While  the  signal  may  be 
able  to  be  salvaged  by  hltering  the  true  gradient  as  well,  this  illustrates  the  effects  of 
the  low  pass  hlter  at  high  velocities  -  signihcant  loss  of  spatial  frequency  information 
in  addition  to  time  delay.  While  no  conventional  vehicles  currently  inhabit  this  flight 
regime,  it  is  important  to  note  the  importance  of  including  the  moving  average  effect 
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that  real  time  signal  generation  will  require.  Additionally,  this  behavior  further  serves 
to  validate  the  overall  modeling  effort.  Addressing  these  phenomena  and  optimizing 
the  LPF  for  real-time  signal  generation  will  challenge  future  nav-grade  gradiometer 
and  map-matching  algorithm  designers. 


Figure  47:  Signal  Comparison:  Excellent  (GGIl)  vs.  Truth,  300m/s,  2500m,  Rough 
Terrain. 

Figure  47  presents  a  closer  look  at  the  behavior  of  an  excellent  signal  when 
compared  to  the  true  gradient.  This  signal  was  generated  by  the  lower  noise  GGI  at 
2500m  HAAT,  a  velocity  of  300m/s  and  over  rough  terrain  (run  28).  Note  that  while 
the  time  delay  is  still  evident,  the  signal  clearly  outlines  the  overall  gradient  contour. 
Also  note  that  the  contour  is  unique  yet  is  low  enough  in  frequency  content  such  that 
the  LPF  does  not  discard  useful  information.  This  gives  foresight  into  the  idea  that 
a  contour  matching  algorithm  may  be  a  good  starting  point  for  design  of  a  gravity 
gradient  map-matching  scheme. 

Figure  48  shows  a  useful  signal  produced  by  the  lower  noise  GGI  in  comparison 
to  the  true  gradient.  This  signal  was  generated  at  10000m  HAAT,  a  velocity  of 


Figure  48:  Signal  Comparison:  Useful  (GGIl)  vs.  Truth,  600m/s,  10000m,  Smooth 
Terrain. 

600m/s  and  over  smooth  terrain  (run  83).  Note  the  overall  flattening  of  the  true 
gradient  when  compared  to  the  previous  case  -  this  serves  to  reduce  the  uniqueness  of 
the  signal.  As  such,  GGIl  produces  a  useful  signal  by  doing  a  fair  job  of  highlighting 
the  true  gradient  contour.  While  changes  changes  aren’t  captured  every  second,  the 
gradient  contour  is  dehned  by  the  signal. 

Figure  49  shows  a  marginally  useful  and  unusable  signal  in  comparison  to  the 
true  gradient.  These  signals  were  generated  at  10000m  HAAT,  a  velocity  of  300m/s 
and  over  smooth  terrain  (runs  82  and  94).  Note  that  the  true  gradient  has  become 
relatively  flat  and  doesn’t  change  more  than  O.SUo  over  the  entire  run.  This  is  a 
difficult  case  for  either  gradiometer.  The  lower  noise  gradiometer  (GGIl)  produces  a 
marginally  useful  signal  by  doing  a  fair  job  of  highlighting  the  true  gradient  contour. 
Even  so,  since  the  contour  doesn’t  change  very  rapidly,  it  is  more  difficult  to  ascertain 
the  exact  shape  of  the  contour.  The  noisier  gradiometer  (GGI2),  however,  produces 
a  signal  that  gives  no  useful  information  about  the  contour  shape.  For  scenarios  with 
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Figure  49:  Signal  Comparison:  Marginally  Useful(GGIl)  vs.  Unusable(GGI2)  vs. 

Truth,  300m/s,  10000m,  Smooth  Terrain. 

relatively  sparse  terrain  (and  relatively  constant  subterranean  features),  it  is  likely 
that  the  gradiometer  must  produce  noise  an  order  of  magnitude  lower  than  GGIl 
(approx.  Q.QlEo  of  error)  for  the  signal  to  be  excellent.  This  is  consistent  with 
results  from  Richeson’s  high  altitude,  hypersonic  cases  [51]. 

This  concludes  the  qualitative  gradiometer  signal  analysis  which  has  shown  that, 
as  generally  expected,  higher  altitudes,  smooth  terrain,  slower  velocities,  and  increased 
noise  can  signihcantly  reduce  the  usefulness  of  the  GGI  signal.  Additionally,  the 
analysis  has  shown  that  LPF  settings  must  considered  and  carefully  chosen  so  as  not  to 
exclude  useful  high  frequency  gradient  information  when  traveling  at  low  altitudes  and 
high  velocities.  For  a  constant  noise  level,  the  largest  contributor  to  signal  degradation 
is  the  overflight  of  relatively  flat  terrain  with  constant  subterranean  density.  The  next 
largest  contributor  to  signal  loss  is  an  increase  in  altitude,  though  slower  velocities 
are  nearly  as  detrimental.  Unfortunately,  the  most  difficult  regimes  to  obtain  a  useful 
signal  are  also  the  regimes  where  most  conventional  military  aircraft  operate  (50- 


300m/s,  5-10A;m).  In  order  to  ensure  an  excellent  to  usable  signal  is  present  for 
all  tested  flight  conditions,  a  gradiometer  with  O.OlEo  of  RMS  error  will  likely  be 
required. 

Navigation  Feasibility  via  Signal  Time  Rate  of  Change  Metries.  Figure  50 
shows  the  signal  usefulness  envelope  based  on  metrics  obtained  from  Richeson’s  sim¬ 
ulations  and  is  a  summary  of  the  results  shown  in  Figures  C.11-C.15  located  in  Ap¬ 
pendix  C.  Recall  that  these  metrics  only  apply  to  GGIl  and  that  the  amount  of  noise 
in  GGI2  gave  Richeson  no  improvement  in  navigation  performance  over  an  unaided 
INS.  Immediately  evident  is  the  overall  signal  degradation  when  compared  to  the  pre¬ 
vious  method  (see  Figure  44),  particularly  in  areas  with  smooth  terrain.  As  before, 
low  terrain  variance  and  increases  in  altitude  are  the  main  contributors  to  signal  degra¬ 
dation.  At  1000m  over  both  terrain  types,  all  velocities  give  signals  that  are  at  least 
marginally  useful,  with  most  being  useful  to  excellent.  However,  the  results  quickly 
change  as  altitude  increase.  The  signal  at  2500m  has  already  drastically  dropped  for 
the  slower  velocity  cases.  For  the  rough  terrain,  the  signal  is  generally  useful,  whereas 
for  the  smooth  terrain,  the  signal  is  only  marginally  useful  for  velocities  greater  than 
300m/s.  At  5000m,  the  signal  has  degraded  signihcantly  -  only  a  marginally  useful 
to  useful  signal  is  obtained  for  the  rough  terrain  case  and  a  largely  unusable  signal 
exists  for  smooth  terrain  overflight.  Note  that  at  5000m,  faster  velocities  (600m/s-|-) 
are  critical  in  order  to  maintain  signal  usefulness  over  rough  terrain  due  to  the  drop  in 
gradient  strength.  The  signal  continues  to  degrade  as  the  gradients  further  attenuate 
at  10000  and  20000m.  The  rough  terrain  signal  becomes  marginally  useful  while  the 
smooth  terrain  signal  is  completely  unusable.  If  it  is  again  assumed  that  most  military 
aircraft  will  be  operating  in  the  50-300m/s  and  5-10A;m  region  of  the  envelope,  the 
sensed  signal  is  marginally  useful  (rough  terrain)  to  unusable  (smooth  terrain). 


Altitude  (km)  Altitude  (km) 


Unusable 


Figure  50:  GGIl  Signal  Glassification  Summary  -  Signal  Time  Rate  of  Ghange 
Method. 
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Figure  51:  Signal  Comparison:  Marginally  Useful  (GGIl)  vs.  Unusable  (GGI2)  vs. 

Truth,  300m/s,  10000m,  Rough  Terrain,  Signal  Rate  of  Ghange  Metrics. 

For  illustrative  purposes.  Figure  51  shows  the  signal  required  for  the  heart  of 
the  “marginally  useful”  envelope  (300m/s  ground  velocity,  10000m  altitude,  rough 
terrain).  This  marginally  useful  signal,  produced  by  GGIl,  is  shown  in  red.  For 
comparative  purposes,  the  signal  from  the  noisier  gradiometer  (GGI2)  is  also  shown. 
Note  how  close  the  GGIl  signal  appears  to  the  noise  free  signal,  yet  it  still  produces  a 
“marginally  useful”  signal.  Also  note  that  GGI2  is  able  to  capture  the  general  shape 
of  the  contour,  albeit  not  as  well  as  GGIl.  This  gives  rise  to  the  question  of  whether 
a  revised,  gradiometer  specific  contour  based  map-matching  algorithm  could  provide 
substantial  benefits  in  navigation  accuracy. 

Figure  52  shows  the  signal  from  the  gradiometers  in  the  heart  of  the  unusable 
portion  of  the  envelope  (300m/s  ground  velocity,  10000m  altitude,  smooth  terrain). 
Recall  that  GGIl  produced  a  marginally  useful  signal  in  this  same  scenario  when  qual¬ 
itative  metrics  were  used  (see  Figure  49).  Evident  is  the  fact  that  the  signal  changes 
very  little  in  time.  While  slight  contour  changes  may  be  measured  by  GGIl,  those 
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Figure  52:  Signal  Comparison:  Unusable  (GGIl  and  2)  vs.  Truth,  300m/s,  10000m, 
Smooth  Terrain,  Signal  Rate  of  Ghange  Metrics. 

changes  are  likely  to  take  several  minutes  to  appear.  Gradiometer  accuracy  would 
need  to  improve  at  least  one  order  of  magnitude  over  that  of  GGIl  (approximately 
O.OlUo  error)  before  a  usable  navigation  signal  could  be  obtained  from  this  worst-case 
scenario.  While  these  results  differ  slightly  from  the  qualitative  study,  results  from 
both  cases  suggest  that  a  O.OlUo  gradiometer  will  be  required  to  ensure  a  useful  signal 
for  all  portions  of  the  envelope,  particularly  for  the  conditions  that  most  conventional 
military  aircraft  operate  in.  That  stated,  a  gradiometer  with  performance  similar  to 
GGIl  could  be  used  for  an  airborne  navigation  feasibility  flight  test  or  demonstration 
if  flown  at  relatively  low  altitudes  and  over  rough  terrain. 

To  test  the  recommended  gradiometer  specihcation  of  O.OlUo,  a  run  was  flown 
at  20000m  over  smooth  terrain  at  a  velocity  of  50m/s.  The  results,  shown  in  Figure 
53,  show  the  signal  produced  by  this  gradiometer  at  the  worst  possible  test  conditions. 
The  plot  shows  that  signal  is  marginally  useful  as  it  is  able  to  highlight  the  overall 
gravity  gradient  contour,  but  it  does  take  time  to  do  so.  Also  note  the  gradient  scale 
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in  the  figure.  Not  only  does  the  GGI  have  to  provide  a  signal  with  O.OlEo  of  error, 
but  the  gradient  maps  to  which  this  signal  will  be  matched  must  contain  considerably 
less  error.  These  are  the  challenges  that  must  be  met  to  ensure  a  useful  signal  in  all 
flight  conditions. 


Figure  53:  Signal  Gomparison:  Ultra  Low  Noise  GGI  (O.OlEo)  vs.  Truth,  Worst  Gase 
Scenario  -  50m/s,  20000m,  Smooth  Terrain. 


The  Tanker  Effect.  Figure  54  captures  the  overall  effect  of  a  large  tanker 
aircraft  in  the  presence  of  the  GGI-carrying  aircraft  flying  at  150m/s  and  5000m 
HAAT  (runs  121-122).  Only  one  figure  is  presented  because  the  trends  are  the  same 
for  each  run  -  the  tanker  essentially  biases  the  signal  by  approximately  0.45Eo  and 
adds  a  slight  amount  of  white  noise  (ct=0.18Eo)  corresponding  to  relative  motion 
between  the  tanker  and  GGI-aircraft  as  the  formation  is  maintained. 

In  no  case  does  the  tanker  change  the  usefulness  of  the  signal  produced  by  the 
two  GGIs.  Should  GGI  performance  improve  an  order  of  magnitude  over  GGIl  (i.e.  a 
O.OlEo  GGI),  the  presence  of  the  tanker  may  begin  to  hamper  map-matching  perfor- 
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Figure  54:  Tanker  Effect,  GGIl,  5000m,  150m/s. 
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mance,  particularly  over  smooth  terrain.  That  stated,  a  GGI  specihc  map-matching 
algorithm  should  be  inherently  designed  to  reject  biases.  It  could  be  argued  that 
the  presence  of  the  tanker  is  akin  to  a  self-gradient  and  could  be  removed  similarly. 
However,  without  broadcasting  an  active  signal,  the  exact  position  of  the  tanker  is  un¬ 
known.  Plots  of  the  remaining  tanker  test  cases  (runs  123-128)  are  shown  in  Figures 
G.16-G.18  located  in  Appendix  G. 

Aircraft  Terrain  Avoidance  Feasibility  Results 

The  results  for  the  hrst  terrain  avoidance  method  will  now  be  presented.  Recall 
that  no  prior  positional  information  is  known  and  a  GGI  with  a  IHz  gradient  produc¬ 
tion  rate  and  0.2Hz  cutoff  is  being  used  as  the  sole  device  in  an  attempt  to  provide 
a  consistent  terrain  avoidance  warning.  It  should  be  noted  that  much  of  the  signal 
for  the  relatively  small  obstacles  used  in  these  simulations  is  in  the  higher  frequency 
(shorter  wavelength)  portion  of  the  spectrum.  As  such,  the  cutoff  frequency  of  the 
LPF  may  present  an  issue  when  trying  to  detect  the  smaller  obstacles.  The  aircraft  is 
flying  at  50m/ s  and  is  level  on  a  plane  located  10m  below  to  the  top  of  the  obstacle 
(see  Figure  35).  Both  the  T^x  signal  and  signal  time  rate  of  change  are  examined  in 
an  attempt  to  determine  a  usable  threshold  that  indicates  imminent  terrain  collision. 

Figure  55  summarizes  the  GGI  signal  and  signal  time  rate  of  change  for  the 
5  runs  against  obstacles  of  varying  size.  Results  from  the  individual  runs  can  be 
found  in  Appendix  G,  Figures  G.19-G.23.  The  analysis  begins  by  examination  of 
the  simulation  using  the  smaller  obstacles  (25  and  50m).  Given  the  requirement  for 
a  warning  1.5s  prior  to  impact  and  a  gradient  production  rate  of  IHz,  the  worst 
possible  time  for  an  update  is  approximately  2.49s  prior  to  the  impact.  An  update 
2.49s  prior  to  impact  maximizes  the  time  spent  without  an  update  by  putting  the 
next  update  at  1.49s  until  impact  -  too  late  given  the  threshold  used  for  this  study.  In 
other  words,  the  update  at  2.49s  prior  to  impact  must  contain  enough  information  to 
trip  a  warning  to  the  operator.  The  figure  clearly  illustrates  that  the  hltered  signal 
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Time  to  Terrain  impact  (s) 

(a)  GGI  Signal. 


Figure  55:  Terrain  Avoidance  Scenario  Summary,  GGI  signals,  50  m/s.  No  Noise. 
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lags  the  true  gradient  considerably  and  that  since  most  of  the  smaller  obstacles’  signal 
intensity  is  in  the  high  freqnency  area  of  the  spectrum,  it  is  hltered  out.  Thus,  there 
is  no  usefnl  information  abont  the  obstacles  at  2.49s.  Examination  of  the  results  for 
the  simulation  using  a  100m  cubic  obstacle  shows  that  while  there  is  a  slight  rise  in 
the  GGI  signal  at  2.49s,  it  is  likely  to  be  bnried  in  signal  noise.  The  resnlts  for  the 
simnlation  using  a  250m  cnbic  obstacle  show  snbstantial  improvements  in  gradient 
strength  and  time  of  detection  over  previous  runs.  Finally,  as  expected,  the  results 
for  the  simulation  using  a  500m  cnbic  obstacle  show  even  larger  valnes  and  rates  of 
change  of  the  sensed  gradient.  Unfortnnately,  the  plots  also  clearly  show  that  there  is 
no  red  flag  signal  or  signal  rate  of  change  that  signihes  imminent  terrain  impact.  The 
lag  and  loss  of  short  wavelength  information  from  the  LPF  are  also  evident.  These 
results  are  not  entirely  nnexpected  dne  to  the  inverse  natnre  of  using  gravity  gradients 
for  what  is,  in  essence,  ranging  information.  There  are  an  inhnite  nnmber  of  possible 
obstacles  with  different  densities  and  locations  that  conld  provide  the  same  signal. 
Also  note  that  these  runs  were  accomplished  at  only  50m/s,  a  relatively  low  velocity 
in  comparison  to  many  terrain  following  aircraft  that  routinely  fly  low-level  rontes 
at  200-|-m/s.  At  this  point,  the  investigation  moves  to  an  examination  of  the  actnal 
gradients  in  order  to  determine,  even  if  a  perfect  gradiometer  existed,  if  there  is  some 
signal  threshold  that  gives  obstacle  ranging  information. 

Fignre  56  summarizes  the  true  T^x  gradient  and  gradient  rate  of  change  prodnced 
by  the  5  obstacles  as  the  simnlation  grid  is  traversed.  The  plots  clearly  show  that, 
even  with  a  perfect  gradiometer  capable  of  measnring  the  true  signal,  there  is  no 
clear  signal  or  signal  rate  of  change  threshold  that  signihes  imminent  terrain  impact 
for  these  simple  scenarios.  While  the  measnred  gradient  time  rate  of  change  contains 
more  information  earlier  in  time  than  the  measnred  gradients,  there  is  no  nniqneness 
based  on  time  to  impact.  If  the  threshold  was  arbitrarily  set  at  some  valne,  the  user 
will  get  false  alarms  for  larger  obstacles  and  no  alarm  at  all  for  small  objects.  Avoiding 
a  big  monntain  does  no  good  if  the  user  subsequently  impacts  a  small  monnd  of  dirt. 
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Rate  of  Change  (Eo/s) 


This  is  the  crux  of  GGI-based  terrain  avoidance  -  without  additional  information, 
ranging  information  cannot  easily  be  determined  with  only  a  gradiometer. 


Time  to  Terrain  Impact  (s) 


(a)  True  Gradient. 


Figure  56:  Terrain  Avoidance  Scenario  summary,  True  Gradients,  50  m/s. 


Recall  that  the  second  method  of  terrain  avoidance  assumes  that  exact  positional 
information  is  known  and  that  gradient  and  terrain  database  maps  are  available. 
The  role  of  the  gradiometer  then  is  to  detect  obstacles  that  were  un-modeled  on  the 
gradient  and  terrain  database  maps.  A  large  water  tower  was  chosen  to  represent 
a  best  case  un-modeled  obstacle.  That  is,  this  object  is  the  largest  likely  to  be 
constructed  quickly  enough  to  avoid  being  included  in  NOTAMs  or  in  intelligence 
reports.  If  the  gradiometer  cannot  warn  the  user  of  the  presence  of  an  obstacle  this 
size,  scenarios  with  smaller  objects  such  as  communication  towers  will  be  even  less 
successful.  To  better  understand  the  signal  structure  of  the  obstacle,  the  analysis  of 
this  terrain  avoidance  method  begins  by  examining  the  T^x  gradient  produced  by  the 
water  tower  terrain  anomaly  in  the  spatial  frequency  domain. 


Figure  57:  Water  Tower  Spectrum 

Figure  57  shows  the  absolute  amplitudes  of  the  gradient  as  a  function  of  spa¬ 
tial  frequency  (for  continuity,  this  spectrum  was  calculated  on  a  horizontal  plane  5m 
above  the  top  of  the  tower).  The  hgure  shows  that  most  of  the  signal  intensity  is 
located  in  the  0. 005-0. 03cj/c/m  spatial  frequency  region  (33-200m  wavelengths).  This 
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immediately  raises  concern  because,  even  at  relatively  slow  speeds  (~  50m/s),  the 
minimum  wavelength  allowed  by  a  LPF  with  a  0.2Hz  cutoff  frequency  is  approxi¬ 
mately  250m  (see  Table  9).  This  mirrors  the  phenomenon  seen  in  the  25  and  50m 
obstacles  from  the  previous  analysis. 


(a)  Signal. 


(b)  Signal  rate  of  change. 


Figure  58:  Water  Tower  Scenario,  Txx  vs  time,  50m/s,  No  Noise. 
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The  results  for  the  water  tower  run  are  shown  in  Figure  58.  As  in  the  25  and 
50m  obstacle  results,  a  gradiometer  providing  gradients  at  IHz  and  hltering  the  signal 
above  0.2Hz  results  in  the  water  tower  gravitational  gradients  never  being  sensed  by 
the  gradiometer  until  after  impact.  Also  shown  are  the  true  gradients  prodnced  by  the 
water  tower.  If  the  assnmption  of  a  perfect  gradiometer  is  made,  the  water  tower  does 
produce  a  signal  of  approximately  l.OEo  at  1.5s  to  tower  impact.  Again  the  question 
arises:  is  the  anomalous  gradient  from  a  small  tower  close-by  or  a  larger  un-modeled 
object  in  the  distance  (i.e.  the  inverse  problem).  Even  if  ranging  information  were 
able  to  be  determined  (perhaps  by  using  the  terrain  avoidance  method  proposed  by 
Jircetano  [46]),  the  question  of  reqnired  gradiometer  performance  arises.  To  determine 
these  valnes,  the  gradiometer  requirements  for  sensing  most  of  the  example  water 
tower’s  spectrum  (down  to  33m  wavelengths)  at  different  velocities  are  listed  in  Table 
14. 


Table  14:  Gradiometer  Reqnirements  to  sense  Water  Tower 


Velocity  (m/s) 

Est.  Gradient  Production  Rate(Hz) 

Cutoff  Frequency  (Hz) 

50 

7.5 

1.5 

100 

15 

3.0 

150 

23 

4.5 

300 

45 

9.0 

600 

91 

18 

1200 

182 

36 

To  further  support  these  estimated  gradiometer  requirements,  it  is  assnmed  that 
additional  information  is  snpplied  via  a  gravimeter  in  a  scheme  similar  to  the  Lockheed 
Martin  UGM.  Under  the  assnmption  that  the  UGM  provides  adeqnate  snbmarine 
terrain  avoidance  capability  and  that  the  gradiometer  used  within  the  UGM  provides 
gradients  at  IHz  with  a  cutoff  frequency  of  0.2Hz  and  is  traveling  at  20knots  (lOm/s), 
gradients  are  produced  every  10m  and  wavelengths  greater  than  50m  are  sensed.  By 
converting  these  metrics  into  equivalent  airborne  platform  reqnirements  based  on 
velocities.  Table  15  snmmarizes  the  eqnivalent  gradiometer  update  rate  and  cntoff 
freqnency  reqnirements. 
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Table  15:  Predicted  Gradiometer  Requirements  for  UGM  Style  Terrain  Avoidance 


Velocity  (m/s) 

Gradient  Production  Rate  (Hz) 

Cutoff  Frequency  (Hz) 

50 

5 

1 

100 

10 

2 

150 

15 

3 

300 

30 

6 

600 

60 

12 

1200 

120 

24 

It  must  be  stressed  that  these  update  rates  and  cutoff  frequency  specihcations 
must  be  met  without  an  increase  in  noise.  That  is,  the  gradiometer  must  have  massive 
improvements  in  gradient  production  rate  and  LPF  cutoff  frequency,  yet  produce  sub- 
Eotvos  noise  levels  in  the  signal.  In  the  300m/s  scenario,  the  necessary  gradiometer 
noise  spectral  density  for  O.lEo  of  noise  would  be  approximately  {i.{i2bEo/\/Hz,  valid 
up  to  half  the  sampling  rate  of  ?)QHz.  This  is  well  beyond  any  estimated  performance 
levels  for  future  airborne  gradiometers.  These  results  highlight  the  signihcant  chal¬ 
lenge  of  using  a  gradiometer  based  system  for  aircraft  terrain  avoidance.  Based  on 
the  results  from  both  methods,  it  has  been  shown  that  gradiometer  based  airborne 
terrain  avoidance  is  unfeasible  in  the  near  future. 

Results  Summary 

It  has  been  shown  that  if  airborne  GGIs  can  approach  O.OlEo  standard  devia¬ 
tion  of  noise,  a  signal  strong  and  unique  enough  for  map-matching  exists  for  all  tested 
flight  conditions  except  for  those  involving  extremely  high  velocities  (1200m/s)  and 
very  low  altitudes  (1000m).  This  is  provided  that  accurate  gravity  gradient  databases 
exist.  The  necessity  of  such  a  gradiometer  is  driven  by  the  fact  that  gradient  unique¬ 
ness  falls  rapidly  as  the  host  vehicle  traverses  smooth  terrain,  increases  altitude  or 
decreases  velocity.  To  adequately  measure  a  gravity  gradient  contour  in  a  worst  case 
scenario,  the  noise  level  of  the  instrument  must  be  an  order  of  magnitude  lower  than 
those  of  gradiometers  projected  to  be  available  within  10  years.  Also,  the  assumption 
of  the  existence  of  accurate,  high  resolution  gradient  maps  cannot  go  unchallenged. 
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While  map  availability  is  rapidly  increasing  through  a  variety  of  survey  methods, 
worldwide  coverage,  particularly  in  remote  or  hostile  areas,  will  likely  remain  a  con¬ 
tingency  for  quite  some  time.  While  it  is  unlikely  that  this  navigation  method  can 
be  widely  employed  using  a  gradiometer  available  within  10  years,  the  signal  analysis 
has  also  shown  that  a  GGI  with  O.lEo  of  noise  could  be  used  in  a  ground  or  flight 
test  demonstration  of  the  technology.  This  assumes  that  relatively  accurate  gravity 
gradient  maps  of  the  area  exist  and  the  test  is  done  at  slow  velocities,  low  altitudes, 
and  over  rough  terrain.  Success  in  such  a  test  could  open  the  door  for  more  research 
and  development  funding. 

The  GGI-based  terrain  avoidance  studies  have  painted  a  relatively  bleak  picture 
for  the  method  as  it  pertains  to  conventional  aircraft.  Due  to  bandwidth,  gradient 
production  rate  limitations  and  the  inverse  nature  of  the  problem,  a  noise  free  gra¬ 
diometer  failed  to  provide  any  useful  terrain  avoidance  information.  These  tests  were 
designed  as  a  best  case  approach  to  solving  the  problem  and  proved  that  no  consistent 
signal  threshold  for  imminent  terrain  impact  exists.  Should  the  inverse  problem  of 
ranging  dangerous  obstacles  be  solved  via  a  novel  method  or  increased  observabil¬ 
ity,  real  world  application  of  such  a  scheme  would  inevitably  demand  extremely  high 
gradiometer  performance  -  the  likes  of  which  have  not  been  mentioned  in  open  litera¬ 
ture.  That  stated,  research  into  methods  of  GGI-based  or  assisted  terrain  avoidance 
should  not  be  abandoned  as  they  could  eventually  be  used  in  applications  such  as 
cave  navigation  or  navigation  through  indoor  environments.  With  the  completion 
of  the  results  and  analysis  of  this  study,  the  focus  now  turns  to  recommending  the 
future  steps  needed  to  make  this  magnihcent  instrument  the  game  changer  it  has  the 
potential  to  be. 
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V.  Discussion 


Research  Summary 

This  research  represents  the  hrst  step  in  a  multi-phase  approach  towards  de¬ 
velopment  of  a  gravity  gradiometer-based  aircraft  navigation  and  terrain  avoidance 
system.  By  generating  gravity  gradient  data  from  terrain  and  underlying  geology, 
a  realistic  representation  of  gradients  produced  by  the  Earth  was  obtained.  Next, 
using  available  gradiometer  specihcations,  a  signal  processing  approach  was  imple¬ 
mented  to  model  the  effects  of  real-time  gradient  measurement  onboard  an  aircraft. 
The  resulting  signal  was  analyzed  via  metrics  developed  to  rate  signal  strength  and 
uniqueness  at  a  variety  of  representative  flight  conditions.  Based  on  this  compari¬ 
son  between  GGI  signals  and  truth,  map-matching  and,  in  turn,  aircraft  navigation 
feasibility  were  demonstrated.  Additionally,  future  gradiometer  performance  require¬ 
ments  necessary  to  produce  a  usable  signal  over  the  entire  tested  flight  envelope  were 
proposed  and  tested.  Finally,  several  best  case  terrain  avoidance  scenarios  were  de¬ 
vised  to  determine  feasibility  of  using  a  GGI  in  such  a  role.  By  using  time  to  terrain 
impact  metrics,  an  attempt  to  hnd  a  threshold  gradiometer  signal  level  was  made 
for  a  variety  of  obstacles.  While  feasibility  for  GGI-based  terrain  avoidance  was  un¬ 
able  to  be  demonstrated,  some  alternative  methods  and  corresponding  gradiometer 
requirements  were  presented. 

Challenges  and  Limitations 

The  two  key  obstacles  which  must  be  overcome  in  order  to  make  gravity  gradient 
map-matching  aircraft  navigation  systems  a  reality  are  the  performance  of  airborne 
gradiometers  and  the  availability  of  accurate  gravity  gradient  maps.  In  this  research, 
it  was  assumed  that  with  continued  interest  from  the  geophysical,  mining,  and  defense 
industries,  gravity  gradiometers  will  eventually  meet  the  proposed  requirements  for 
aircraft  navigation  feasibility.  However,  they  must  meet  these  requirements  -  mak¬ 
ing  accurate  measurements  at  the  level,  yet  still  be  small  and  light  enough 

to  £t  within  limited  space  inside  an  aircraft.  Also,  they  must  be  robust  enough  to 
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survive  the  military  flight  environment  -  including  rapid  maneuvers,  elevated  g-forces 
and  hostile  environmental  conditions  often  encountered  at  forward-deployed  locations. 
Finally,  increased  gradiometer  sensitivity  does  not  come  without  drawbacks.  Account¬ 
ing  for  self  gradients  produced  by  the  host  vehicle  is  already  a  challenging  prospect 
at  the  lEo  sensitivity  level.  For  this  study,  it  was  assumed  that  self  gradients  were 
always  exactly  known.  Should  a  gradiometer  able  to  sense  changes  at  the  O.OlF'o 
level  be  implemented,  precise  monitoring  and  correction  of  self  gradients  will  be  vital 
to  ensure  proper  measurements.  Fuel  slosh,  the  release  of  stores,  control  surface  and 
even  pilot  movement  will  all  present  gradient  changes  that  must  be  accounted  for. 

It  was  also  assumed  that  the  generated  gravity  gradient  maps  represent  truth 
-  that  is,  each  value  was  exactly  correct  in  magnitude  and  position.  Actual  gravity 
gradient  maps  with  enough  resolution  to  be  useful  in  a  map-matching  algorithm  are 
very  limited  in  quantity,  generally  proprietary  and  contain  GPS  level  position  errors. 
Obtaining  accurate  surveys  over  areas  considered  unfriendly  or  hostile  presents  an 
additional  challenge.  However,  as  more  accurate  ground,  air,  and  space-based  surveys 
take  place  and  gradient  calculation  methods  evolve,  the  resolution  and  availability 
(particularly  in  remote  areas)  of  gravity  gradient  maps  will  signihcantly  increase.  In 
summary,  there  are  significant  challenges  and  limitations  to  overcome  before  GGI- 
based  aircraft  navigation  can  be  a  reality.  Given  continued  research  efforts,  these 
issues  are  certainly  solvable.  Though  unlikely  to  be  completely  addressed  within  10 
years,  it  is  anticipated,  based  on  previous  trends,  that  these  issues  will  largely  be 
solved  in  20-30  years  time. 

GGI-based  airborne  terrain  avoidance  is  the  more  difficult  problem  due  to  the 
inverse  nature  of  ranging  hazardous  terrain  via  gravity  gradiometry.  While  meth¬ 
ods  have  been  proposed  to  solve  this  problem  [46,47],  there  is  very  little  research 
published  in  open  literature  to  back-up  these  ideas.  Under  the  assumption  that  the 
inverse  problem  is  solvable,  the  gradiometer  requirements  for  sensing  obstacles  which 
threaten  low  flying  aircraft  are  dramatically  higher  than  any  current  or  proposed  GGI 
specihcations.  That  is  not  to  say  that  GGI-based  (or  assisted)  airborne  terrain  avoid- 
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ance  is  impossible,  but  it  is  a  distinctly  greater  challenge  than  aircraft  navigation  via 
gravity  gradient  map-matching. 

Significant  Contributions  and  Insights 

In  regards  to  gravity  gradient  modeling,  this  research  has  clearly  implemented 
methods  that  were  unable  to  be  found  in  previous  navigation  via  gravity  gradiometry 
works.  First,  a  user  friendly  interface  to  generate  gravity  gradient  data  representing 
realistic  values  produced  by  the  Earth  was  developed.  Next,  it  was  proven  that  if  a  bias 
is  acceptable,  Parker’s  method  can  produce  gravity  gradient  data  much  faster  than 
more  rigorous  methods.  Also,  terrain  and  geology  (high  and  low  frequency)  gravity 
effects  were  accounted  for  using  a  combination  of  modeling  techniques.  Finally,  tanker 
aircraft  and  man-made  ground  obstacle  gravitational  gradient  models  were  developed 
as  well. 

While  relatively  simple  gravity  gradiometer  models  have  been  presented  in  pre¬ 
vious  works,  real  time  signal  processing  effects  such  as  low  pass  hltering  are  rarely 
mentioned.  In  most  previous  studies,  it  was  generally  assumed  that  airborne  gra- 
diometers  provide  truth  measurements  plus  white  noise.  This  is  not  entirely  true 
as  the  assumption  fails  to  account  for  phase  lag  and  time  delay  implications  from 
the  LPF.  In  other  words,  extensively  processed  post-flight  data  was  assnmed  to  be 
available  instantly.  For  a  more  realistic  approach,  this  research  developed  and  vali¬ 
dated  two  airborne  gradiometer  models  reqnired  to  provide  real-time  gravity  gradient 
measnrements  for  use  in  a  map-matching  navigation  algorithm  and  terrain  avoidance 
scenarios.  Also,  the  effects  of  time-delay,  phase  distortion  and  loss  of  high  frequency 
information  caused  by  the  low  pass  hlter  were  examined  and  found  to  likely  degrade 
map-matching  performance  -  especially  at  extremely  low  altitudes  and  high  veloci¬ 
ties.  Finally,  fntnre  navigation-specihc  gradiometer  design  targets  were  proposed  and 
validated. 
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Significant  contributions  to  the  study  of  passive  aircraft  navigation  via  a  GGI- 
based  map-matching  system  are  now  presented.  This  research  developed  two  methods 
of  GGI  signal  usefulness  classihcation  and  demonstrated  that  unique  gravity  gradient 
signals  do  exist  and  can  be  measured  by  an  airborne  gradiometer.  The  study  also 
investigated  the  GGI  signal  in  regimes  previously  not  examined  and  found  that  signal 
levels  are  highly  dependent  on  terrain  uniqueness  and  host-vehicle  altitude  and  ve¬ 
locity.  As  such,  low  terrain  variance,  high  altitude,  and  low  velocities  are  detrimental 
to  the  signal.  Unfortunately,  legacy  military  aircraft  flight  regimes  generally  provide 
the  weakest  signal.  Additionally,  this  study  was  the  hrst  to  examine  the  effects  of 
formation  flight  on  the  GGI  signal.  It  was  shown  that  placing  the  GGI  aircraft  next 
to  a  large  tanker  aircraft  gave  a  bias  and  small  amount  of  noise  to  the  GGI  signal. 
While  it  did  not  signihcantly  hamper  the  signal’s  usefulness  for  the  tested  GGIs,  for¬ 
mation  flight  may  reduce  map-matching  performance  if  a  more  sensitive  gradiometer 
is  used,  particularly  in  areas  with  smooth  terrain.  Also,  it  was  shown  that  future 
GGIs  will  need  to  achieve  error  levels  on  the  order  of  O.OlEo  for  navigation  useful¬ 
ness  over  most  of  the  tested  flight  envelopes  (50-1200m/s  velocity,  l-20km  altitude, 
rough  and  smooth  terrain).  However,  should  a  ground  or  flight  test  demonstration  of 
the  navigation  technology  be  desired,  gradiometers  currently  in  flight  test  for  mining 
industries  can  provide  a  useful  navigation  signal  for  lower  flying  vehicles  traversing 
rough  terrain. 

While  previous  works  involving  GGI-based  terrain  avoidance  are  few  in  number, 
this  research  sheds  some  insight  as  to  why.  After  development  of  a  GGI-based  terrain 
avoidance  test  methodology  by  modeling  a  variety  of  obstacles  which  are  hazardous  to 
low  flying  aircraft,  it  was  demonstrated  that  obstacle  range  calculation  is  an  inverse 
problem  -  no  imminent  terrain  impact  signal  threshold  exists.  It  was  also  shown 
that,  should  the  inverse  problem  be  solved,  major  GGI  improvements  are  needed 
for  terrain  avoidance  feasibility.  For  example,  a  45x  improvement  over  current  GGI 
gradient  production  rate  and  LPF  cutoff  frequency  is  needed  for  an  aircraft  traveling 
300m/s  to  sense  small  obstacles.  Additionally,  the  GGI  must  meet  these  specihcations 
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with  no  additional  noise.  The  key  driver  for  these  steep  gradiometer  requirements  is 
the  fact  that  obstacles  most  dangerous  to  low  flying  aircraft,  such  as  towers  and  water 
tanks,  have  a  high  frequency  signal  easily  hltered  out,  buried  in  noise,  or  masked  in 
the  signal  from  surrounding  terrain. 

Recommendations  for  Future  Research 

With  the  initial  feasibility  investigation  complete,  this  section  highlights  future 
research  recommendations.  While  these  recommendations  are  not  ordered  by  priority 
(all  must  eventually  be  addressed),  the  next  logical  step  for  researchers  is  the  devel¬ 
opment  of  a  gravity  gradiometer  specihc  map-matching  algorithm.  Once  developed, 
the  effort  to  integrate  the  algorithm  into  a  navigation  system  can  begin.  This  is  crit¬ 
ical  as  it  will  ultimately  provide  quantitative  navigation  performance  results  that  can 
be  used  to  validate  and  rehne  the  feasibility  metrics  used  in  this  research.  Future 
research  recommendations  are  as  follows: 

•  Continue  development  of  a  contour  based  gravity  gradient  map-matching  algo¬ 
rithm. 

•  Integrate  map-matching  algorithm  into  an  aircraft  navigation  system  and  de¬ 
termine  qualitative  navigation  performance. 

•  Investigate  the  navigation  performance  of  an  INS  using  a  GGI  for  gravity  com¬ 
pensation  and  position  updates  provided  via  a  map-matching  scheme. 

•  Rehne  gradiometer  model  and  examine  effects  of  different  hlter  types. 

•  Further  validate  and  rehne  the  signal  usefulness  metrics  proposed  in  this  re¬ 
search. 

•  Further  validate  navigation  grade  gradiometer  requirements  proposed  in  this 
research. 

•  Analyze  navigation  performance  when  position  and  measurement  errors  are  em¬ 
bedded  within  gravity  gradient  maps  used  as  truth  sources. 
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•  Generate  gravity  gradient  maps  using  methods  that  include  density  anomalies 
within  the  terrain  and  geology. 

•  Model  low  frequency  gravitational  effects  via  the  higher  resolution  EGM2008 
model  (in  lieu  of  EGM96). 

•  Pursue  methods  to  solve  the  inverse  problem  of  terrain  avoidance. 

Conclusion 

While  this  research  has  shown  that  major  improvements  in  GGI  gradient  pro¬ 
duction  rate  and  bandwidth  are  needed  before  GGI  assisted  terrain  avoidance  can 
be  realistically  considered,  it  has  also  proven  that  the  fundamentals  of  using  a  mod¬ 
ern  gravity  gradiometer  as  the  foundation  of  a  completely  passive,  precision  aircraft 
navigation  system  are  sound  and  that  the  method  is  entirely  feasible.  By  modeling 
sensors  that  provide  real-time  measurement  of  gravity  gradients  and  then  analyzing 
the  resulting  signals  via  several  metrics,  it  was  determined  that  a  specialized  GGI 
can  provide  a  signal  strong  and  unique  enough  for  map-matching  utility.  While  the 
limiting  factors  have  been  mentioned,  they  are  certainly  conquerable  given  more  time. 
With  the  results  of  this  study  and  the  efforts  of  previous  researchers,  along  with  the 
contributions  from  the  gravity  gradiometer  development  community,  the  foundation 
for  a  passive,  essentially  unjammable,  precision  aircraft  navigation  system  has  been 
laid.  Forty  years  ago,  airborne  gravity  gradient  surveys  would  never  work.  Today, 
these  once  impossible  surveys  can  be  flown  daily.  Forty  years  from  now,  gravity 
gradiometers  may  very  well  form  the  backbone  of  aircraft  navigation  systems. 
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Appendix  A.  Matlab®  Code 


Listing  A.l: 

7o  Gravity  Gradient  Signal  and  "Truth"  Calculation  ”/« 


y,  based  on  Parker  methods  and  EGM96  .  "/« 

7.  7. 

y,  Marshall  Rogers  2008  °L 

6  y.  7. 


close  all 
clear  all 
11  clc 

G  =  6.67E-11;  yoUniversal  Gravity  Const. 
p  =  2670;  7oAverage  terrain  density 

16  Eotvos  =  lE-9;  “/.use  to  convert  units  to  Eotvos 

M= input (' Rough  or  Smooth  Terrain  [R/S] ? ' , ' s ' ) ; 

y,  pulling  in  elevation  data,  use/add  different  areas  if  desired 
if  ((M=='R' )  I  (M=='r' )) 

21  lat  =  [35  37];  y.NS  Geodetic  coordinates  (WGS84  reference  ... 

ellipsoid)  of  grid 

long=[-122  -120];  7,EW  Geodetic  coordinates  (WGS84  reference 
ellipsoid)of  grid 

[Z ,  refvec]  =  dted  (  '  CA_elev2  '  ,  l,lat,long);  7,  read  in  rough 
area  dted  data,  Z  is  referenced  to  MSL  (the  geoid) 
elseif  ( (M== ' S ' ) I (M== ' s ' ) ) 

lat=[35  37];  y.NS  Geodetic  coordinates  (WGS84  ... 
reference  ellipsoid)  of  grid 

26  long=[-90  -88];  7oEWGeodetic  coordinates  (WGS84  ... 

reference  ellipsoid)  of  grid 
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[Z  ,  refvec]  =  dted  (  '  AK_elev  '  ,  l,lat,long);  7,  read  in  . 
smooth  dted  data,  Z  is  referenced  to  MSL  (the  geoid 
) 

else 

f printf ( ' Wrong  Answer ! ' ) ; 
break ; 

31  end 

figure  7o  plot  the  terrain  contour  map 

geoshowCZ, refvec , 'Display Type' , 'texturemap') ,  colorbar 
X lab el ( ' Longitude '  ,  'Fontsize'  ,20)  ,  ylabel('Latitude'  ,  'Fontsize'  ... 
,20) 

36  setCgca,  'Font size', 20); 
vcb  =  colorbar ; 

set (get (vcb ,  ' Ylabel ' )  ,  ' String '  ,  ' Terrain  Height  (m)  '  ,  'FontSize'  ,24) 
set (vcb,  'Fontsize'  ,20) 

41  alt = input (' Input  altitude  (meters  height  above  average  terrain)  = 
') ; 

track_start=input (' Track  Start  [lat  long]  =  '); 

track_end=input (' Track  End  [lat  long]  =  '); 

vel= input (' Input  Velocity  (meters/sec)  =  '); 

update_rate=input (' Enter  GGl  Update  Rate  (sec)  =  '); 

46  f ilt er_ cut  of f  =  input (' Enter  GGl  Low  Pass  Filter  Cutoff  Frequency  (. 
Hz)  =  '); 

bandwidth  =  f ilter_cutof f  ; 

NSD= input (' Enter  GGl  Noise  Spectral  Density  (E/(Hz~l/2))  =  '); 
f i lename 1 = input (' Input  Filename  for  GGl  Signals  (use  .mat  ... 
extension):  ','s'); 

f i lenameS = input (' Input  Filename  for  True  Map  (use  .mat  extension): 
','s'); 

51  alt  =mean2  ( Z )  +  alt  ;  7,  Height  above  average  terrain  calc 
alt_EGM96=alt+mean2 (Z) ; 

7o7o  Computer  Earth  '  s  Radius  at  central  grid  point 
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a  =  6378137.0;  %  earth  semimajor  axis  in  meters 

56  f  =  1/298.257223563;  "L  reciprocal  flattening 
e2  =  2*f-f~2;  "L  eccentricity  squared 

lat_middle=(min(lat)+max (lat) ) /2; 

lat  _middle2  =  at  and  (  (  1  -  e2  )  *  t  and  ( lat  _middle  )  )  ;  "/convert  geodetic  lat 
to  geocentric  lat 

61  long_middle  =  ( min ( long )+max(long))/2; 

N=a/ sqrt (l-e2*sin(lat_middle2) "2) ; 

X_ECEF=N*cosd(lat_middle2)*cos (long_middle) ; 

Y_ECEF=N*cos (lat_middle2) *sin (long_middle) ; 

Z_ECEF=(N*(l-e2) ) *sin(lat_middle2) ; 

66  Rl  =  sqrt  (X_ECEF~2  +  Y_ECEF'2+Z_ECEF'2)  ;  "/.  Earth's  radius  at  grid  ... 
midpoint 

dlong  =  max  (  abs  (  long  )) -min  (  abs  (  long  ))  ;  "/distance  longitude 
dlat  =max  (  lat  ) -min  (  lat  )  ;  "/  distance  latitude 

x_dist=Rl*(pi/180) *dlong*cosd(lat_middle) ; 

71  y _di St =R1 * ( pi / 180 ) * dlat ; 

y_  int  =  R1  *(  pi  / 180  )*(  3/3600 )  ;  °/3  arc  sec  spacing  (DIED  Level  1)  ... 

3/3600 

x_int=cosd(lat_middle)*y_int ; 

Y1=[0: y_int : y_dist] ; 

X1=[0: x_int : x_dist] ; 

76  [X2 ,Y2]=meshgrid(Xl , Yl) ; 

x_track=-Rl*(pi/180) *(abs (track_end (2) ) -abs (track_start (2) ) ) *cosd( 
lat_middle ) ; 

y_track=Rl*(pi/180) *(abs (track_end (1) ) -abs (track_start (1) ) ) ; 
track_ angle  =  (180/pi)*atan2(y_track ,x_track)  ; 

81 

"/  setup  initial  position  and  velocity  for  simulink 
x_vel= vel * cosd ( tr ack_angle ) ; 
y_vel= vel * s ind ( tr ack_angle ) ; 
vi=[x_vel  y_vel  0]; 
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86  speed=norm ( vi ) ; 

sim_time=floor(sqrt(x_track~2+y_track~2)/vel) ; 

ipos=[x_dist-Rl*(pi/180)*(abs(track_start(2))-min(abs(long)))*cosd 
(lat_middle)  y_dist -R1 * ( pi / 180 ) * abs (max ( lat ) - 1 r ack_st art  (  1 ) )  .. 

0]  ;  "/o  initial  position  based  on  input 

figure  7o  plot  the  track  onto  the  terrain  contour  map 
91  geoshow (Z , ref vec DisplayType texturemap ') ,  colorbar 

xlabelC'Longitude'  ,  'Fontsize'  ,20)  ,  ylabelC'Latitude'  ,  'Fontsize'  ... 
,20) 

setCgca,  'Fontsize'  ,20)  ; 
vcb  =  colorbar ; 

set (get (vcb ,  ' Ylabel ' )  ,  ' String '  ,  ' Terrain  Height  (m)  '  ,  'FontSize'  ,24) 
96  set(vcb,  'Font Size', 20) 
hold  on 

plot ([track_start (2)  track_end(2)] , [track_start (1)  track_end(l)] , ' 
k' , 'Linewidth' ,4) 

“/o"/.  Gradients  by  Parker's  Methods  (Terrain) 

101  del_xl  =  x_int  ;  7,  x  interval 

del_x2  =  y_int  ;  7o  y  interval 
ml=length (XI ) ; 
m2= length ( Y1 ) ; 
pl=[0:l:ml-l]  ; 

106  p2=  [0  :  1  :  m2-l]  ; 

for  ctr  =  1 : ml 

if  pi ( ctr ) < (ml /2 ) -1 

f  l_p  (  ctr  )  =pl  (  ctr  )  /  (  del_xl  *ml  )  ;  7.  spatial  frequency  x 

111  f  2_p  (  ctr  )  =p2  (  ctr  )  /  (  del_x2  *m2  )  ;  7.  spatial  frequency  y 

else 

f  l_p  (  ctr  )  =  (  pi  (  ctr  ) -ml  )  /  (  del_xl  *ml  )  ;  7o  spatial  frequency  x 
f  2_p  (  ctr  )  =  (  p2  (  ctr  ) -m2  )  /  (  del_x2  *m2  )  ;  7o  spatial  frequency  y 

end 

116  end 
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[f  lm_p  ,  f  2m_p  ]  =meshgr  id  (  f  l_p  ,  f  2_p  )  ;  "L  gridded  spatial  freq  data 
f_p  =  sqrt (flm_p . ~2+f2m_p  .  ~2)  ; 


121  sig=0; 

for  ctr=l : 20 

sig=sig+((l./factorial(ctr)) .*(2*pi.*f_p) ."(ctr-2) 
ctr)));  7o  perform  FFT  (see  Jekeli  &  Zhu) 

end 

126  sig(l)=lE15;  “/oprevent  infinite  value 
mu_xx = - ( ( 2* pi ) “2) .*flm_p.~2; 
mu_xy  =  -((2*pi)~2)  .*flm_p.*f2m_p; 
mu_xz=i * ( ( 2* pi ) ~ 2) . * f lm_p . * f _p ; 
mu_yy = - ( ( 2* pi ) “2) .*f2m_p.~2; 

131  mu_yz=i * ( (2* pi ) ~ 2) . * f 2m_p . * f _p ; 
mu_zz = ( ( 2* pi ) "2) .*f_p.~2; 

Txx_parker=(2*pi*p*G . *ifft2 ( mu_xx .*exp(-2*pi*alt .*f_p) 
Txy_parker=(2*pi*p*G . *ifft2 ( mu_xy .*exp(-2*pi*alt .*f_p) 
136  Txz_parker=(2*pi*p*G.*ifft2( mu_xz .*exp(-2*pi*alt .*f_p) 
Tyy_parker=(2*pi*p*G . *ifft2 ( mu_yy .*exp(-2*pi*alt .*f_p) 
Tyz_parker= (2*pi*p*G . * ifft2 ( mu_yz .*exp(-2*pi*alt .*f_p) 
Tzz_parker=(2*pi*p*G . *ifft2 ( mu_zz .*exp(-2*pi*alt .*f_p) 

141  Txx_parker =real ( Txx_parker ) ./Eotvos; 

Txy_parker=real (Txy_parker) . /Eotvos ; 

Txz_parker=real (Txz_parker) . /Eotvos ; 

Tyy_parker=real (Tyy_parker) . /Eotvos ; 

Tyz_parker=real (Tyz_parker) . /Eotvos ; 

146  Tzz_parker =real ( Tzz_parker ) ./Eotvos; 

y,7o  Gradients  from  EGM96  (Long  wavelength  subterranean 
7o  Taken  from  Kiamehr  &  Eshagh  and  Modified 
7o  Could  also  use  geopot97  .  vO  .  4e  .  f  code 


.  *fft2  (  (Z . ' 


. *  sig)  ) 
. *  sig)  ) 
. *  sig)  ) 
. *  sig) ) 
. *  sig)  ) 
. *  sig)  ) 


effects) 
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151  phi_south=min ( lat ) ; 
phi_north=max ( lat ) ; 
lambda_we St =min ( long ) ; 
lambda_east =max ( long ) ; 

phi_step  =30/60 ;  7oEGM96  provides  a  30  arcmin  resolution 
156  lambda_step  =30/60 ;  7oEGM96  provides  a  30  arcmin  resolution 

X_EGM96  =  lambda_west  :  3/3600  :  lambda_east  ;  /osetup  3  arcsec  array  for 
griddata  function 

Y_EGM96  =  phi_south  :  3/3600 ;  phi_north  ;  "/setup  3  arcsec  array  for  ... 
griddata  function 

[X2_EGM96 , Y2_EGM96] =meshgrid (X_EGM96 ,Y_EGM96) ; 

161  f ilename = ' EGM96Gr adient s ' ; 

[Nmax  ,  Ae  ,  GM  ,  C  ,  S  ,  dC  ,  dS]  =Modelread  (' egm_coef  .  ascii ')  ;  7  Read  ... 
spherical  harmonic  model 

[a  ,  b  ,  c  ,  d , g , h , beta , psi , mu , eta]  =  coef f i c ient s ( Nmax +3)  ;  7  calculate  . 

Legendre  coeffs. 

CN=Normal(GM,Ae,Nmax) ; 

166  C (3;  11 , 1) =C (3 : 11 , 1) -CN  (3  ;  11)  '  ;  7  Generation  of  the  Potetial  ... 
Anomaly 

f  id  =  f  open  ( f  ilename  w ')  ;  "7  Opening  a  file  for  the  EGM96  Gradients 

for  phi =phi_ south : phi_step : phi_north 
phigeodetic=phi ; 

171  phi =phi *pi / 180 ; 

"/„  Compute  the  Geocentric  latitude  via  geodetic  latitude 
e2=  .  00669437999013  ;  "/Ist  eccentricity  squared 
phi  =  atan((l-e2)*tan( phi ) )  ; 

176 

°/o  Compute  the  Associated  Legendre  f unctions 
[pnm , dP] =Pnm (phi*180/pi,Nmax+3,Nmax+3) ; 

for  lambda=lambda_west : lambda_step : lambda_east 
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181 


lambda= lambda *pi / 180 ; 


186 


191 


196 


201 


206 


211 


sum  =  0;  7o  Initialize  summations 

suml =0 ; 

sum2  =0 ; 

sum3  =0 ; 

sum4=0 ; 

sum5  =0 ; 

sumN  =  0 ; 

sumdg=0 ; 

sumeta=0 ; 

sumpsi =0 ; 

7,  Comput  at  ion  of  geocenric  distance 
N=Ae/ sqrt (l-e2*sin(phi) ~2) ; 

X_ECEF  =  (N+alt_EGM96) *cos ( phi ) *  cos ( lambda )  ; 
Y_ECEF  =  (N+alt_EGM96) *cos ( phi ) *  sin ( lambda )  ; 
Z_ECEF=( (N+alt_EGM96) *(l-e2) ) *sin(phi) ; 
r=sqrt ( X_ECEF '2+ Y_ECEF ~ 2+ Z_ECEF " 2) ; 


for  n=3 ; Nmax+1 
for  m= 1 : n 


CS  =  (C(n,m)*cos ((m-1) *  lambda ) +S(n,m)*sin((m-1) *. 

lambda ) ) ; 

AA= ( Ae /r ) " (n+2) ; 

AAl=(Ae/r) ~n; 

CSl=(-S(n,m)*cos((m-l)* lambda ) +C(n,m)*sin((m-1) 
♦lambda) ) ; 

PNM=pnm (n , m) ; 

if  ( abs  (m)  -2)  <  0 

PP=(-1) ~(abs(m-4) -l)*pnm(n,abs((m) -4)) ; 
PP1=(-1) "(abs (m-4) -1) ♦pnm (n-l,abs((m)-4)) ; 
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else 


216  PP  =  pnni  (n  ,  abs  (m) -2)  ; 

PPl  =piim  (n-l,abs(m)-2); 

end 

221  if  (abs(m)-l)  <  0 

QQ=(-1) ~ ( abs (m-3) -1) *pnm (n , abs (abs (m) -3) ) ; 
QQ1=(-1) "(abs (m-3) -1) *pnm (n-l,abs(abs(m)... 
-3))  ; 

226  else 

QQ  =  pnm (n, abs (m) -1)  ; 

Q Q 1 =  pnm (n-l,abs(m)-l)  ; 

end 

231  7o  Computing  the  Txx  summation 

suml  =  suml +AA*CS*(a(n , abs (m) ) *PP+b(n, abs (m) ) *pnm 
(n,abs(m))+. . . 

c(n,abs(m)) *pnm (n,abs(m)+2)) ; 

7o  Computing  the  Txy  summation 

236  sum3  =  sum3  + AA*  CS 1 *(d(n,m)*PPl+g(n,m) *pnm (n- 1 , (m) 

) +h (n , m) *pnm (n-l,(m)+2)); 

7o  Computing  the  Txz  summation 

sum4=sum4+AA*CS* (beta(n,m)*QQ+psi(n,m)*pnm(n, (m 
)+l)) ; 

241  7o  Computing  the  Tyz  summation 

sum5  =  sum5  +  AA* CS 1  * (mu (n,m) *QQl  +  eta(n,m) *pnm (n.  .  . 
-1 , (m) +1) )  ; 
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246 


251 


256 


7o  Computing  the  Tzz  summation 
sum2=sum2+ (n* (n+1) ) *AA*CS*PNM ; 

end  7o  of  m 
end  7o  of  n 

7oThe  gravity  gradient  tensor  component  s 

Txx  =  -GM/Ae~3*  suml /Eotvos ; 

Tzz=  GM/ Ae ~3* sum2 / Eotvos ; 

Txy=  (-GM/Ae "3* sum3/Eotvos ) /lO ; 

Tyz=  -GM/Ae ~3* sum4/Eotvos ; 

Txz=  GM/ Ae ~3* sum5 / Eotvos ; 

f  printf  (f  id  ,  '  7og  7og  7.e  7oe  7.e  7oe  7.e  7oe  \n  '  ,  phigeodetic  ,  .  .  . 
lambda*180/pi , . . . 

Txx,-(Txx+Tzz) ,Tzz,Txy,Txz,Tyz) ; 


261  end  7.  of  lambda 

end  7o  of  phi 


fclose(fid)  ; 

U=load ( f ilename ) ; 

266 

7o  interpolate  EGM96  gradients  from  30arcmin  res  to  3arcsec  res 
Txx_EGM96=griddata (U( ; ,2) ,U(: ,1) ,U(; ,3) , X2_EGM96 , Y2_EGM96 , ' v4 ' ) 
Tyy_EGM96=griddata (U (: ,2) ,U(; ,1) ,U(; ,4) , X2_EGM96 , Y2_EGM96 , ' v4 ' ) 
Tzz_EGM96=griddata (U( : ,2) ,U(: ,1) ,U(; ,5) , X2_EGM96 , Y2_EGM96 , ' v4 ' ) 
271  Txy_EGM96=griddata (U( : ,2) ,U(: ,1) ,U(; ,6) , X2_EGM96 , Y2_EGM96 , ' v4 ' ) 
Txz_EGM96=griddata (U (: ,2) ,U(; ,1) ,U(; ,7) , X2_EGM96 , Y2_EGM96 , ' v4 ' ) 
Tyz_EGM96=griddata (U( ; ,2) ,U(: ,1) ,U(; ,8) , X2_EGM96 , Y2_EGM96 , ' v4 ' ) 

7o7o  The  Gradients  ! 

276  Txx=Txx_parker +Txx_EGM96 ; 

Txy=Txy_parker +Txy_EGM96 ; 
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Txz=Txz_parker +Txz_EGM96 ; 

Tyy=Tyy_parker +Tyy_EGM96 ; 

Tyz=Tyz_parker +Tyz_EGM96 ; 

281  Tzz=Tzz_parker +Tzz_EGM96 ; 

7,°/o  Run  Simulink  Model 

simC'GGI',  0 :  s  im_t  ime )  ;  ”/«  used  120s  in  thesis  runs 
movefile ( ' filename2 .mat ' .filenamel) 

286  movefile ( ' filename4 . mat ' , filenames) 


Listing  A. 2: 

7,  This  function  reads  the  spherical  harmonic  model 

3  function  [Nmax , Ae , GM , C , S , dC , dS] =Modelread ( f ilename ) 

fid  =  fopen(filename  ,  'r')  ; 

Al  =  f  scanf  (f  id  ,  '  7.g  7g  7.g  \n',6); 

8 

Nmax  =  Al ( 1)  ; 

Ae=Al (2) ; 

GM  =  A1  (3)  ; 

13  while  (-if  eof  (  f  id  )  ) 

B  =  f  scanf  (fid  7.d  7.d  7.g  7.g  7.g  7og  \n',6); 
n=B(l) ;m=B(2) ; 

C(n+l,m+l)=B(3) ; S (n+1 , m+1) =B (4) ; 

18  dC (n  +  1 ,m  +  l) =B  (5)  ; dS (n  +  1 ,m  +  l) =B (6)  ; 

end 

fclose(fid)  ; 
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Listing  A. 3: 


7.  1 

7,  This  f  unct  ion  computes  the  coefficients  of  the  Legendre  7o 

4  7o  functions  7o 

I  1 

7.  INPUT  7. 

7»  the  maximum  desired  degree  of  geopotential  model  to  be  used7o 

7o  plus  2  -it  is  suggested  to  introduce  higher  values  than  the7o 

9  7o  maximum  degree  of  the  model  .  7o 

7.  7. 

7.  OUTPUT  7. 

7.  7. 

7«  all  of  the  coefficients  of  Legendre  f  unct  ions  needed  for  7o 

14  7o  computing  the  gravity  gradients  7o 

I  1 

7.  REFERENCE  7. 

7o  Petrovskaya,  M.S.  and  A.N.  Vershkov  (2006),  Non/Singular  7o 

7«  expressions  for  the  gravity  gradients  in  the  local  7o 

19  7o  north  -  or  iented  and  orbital  referencse  frames.  Journal  of  7o 

7.  Geodesy,  Vol  80,  117-1277.  7. 

7.  7. 

7.  7. 

7.  by  7. 

24  7.  Mehdi  Eshagh  and  Ramin  Kiamehr  2006  7. 

7.  DivisionofGeodesy  7. 

7.  Royal  Institute  of  Technology  7. 

7.  Stockholm  ,  Sweden  7. 

7.  Email  :  eshaghOkth  .  se  7. 

29  7.  7. 

7.  Modified  by  Marshall  Rogers  2008  7. 

7.  7. 

34 
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function  [a,b  ,  c  ,d,g,li,beta  ,psi  ,mu  ,  eta]  =  coefficients  (N) 
f or  n  =  l ; N 


for  m=l : n 

if  ( ( abs (m-1) ==0) I abs (m-1) ==1) 

39 

a(n , abs (m) ) =70 ; 

b(n, abs (m) ) =(n-l+abs (m-1) +1) *(n-l+abs (m-1) +2) / (2* abs (m 
-1)+1)  ; 

44  c  (n  ,  abs  (m)  )  =  sqrt  (  1  +  A  (  abs  (m- 1 )  ,  0)  )  *  sqrt  (  (n  - 1 )  ~2  -  (  abs  (m  . 

-1) +1) '2) *  .  .  . 


sqrt ( (n-1) -abs (m-1) ) *sqrt (n-l+abs (m- 

-1) +2) /4; 

elseif  (2  <  abs(m-l)  <  n-1) 

a(n, abs (m) )=sqrt (1+A(abs (m-1) ,2) ) *sqrt ((n-1) ~2-(abs (m 
-1) -1) '2) *  .  .  . 


49  sqrt((n-l)+abs(m-l))*sqrt(n-l-abs(m- 

-1)  +2)  /4; 

b(n,abs(m))=((n-l)~2+(m-l) ~ 2+3* (n-1) +2) /2; 

c (n, abs (m) )=sqrt ( (n-1) "2-(abs (m-1) +1) ~2) *sqrt ((n-1) -. . 
abs(m-l))*.  .  . 

54  sqrt ( (n- 1 ) tabs (m- 1 ) +2) /4 ; 

d(n,m)=-(m-l)/4/abs(m-l)*sqrt((2*(n-l)+l)/(2*(n-l)-l)) 


*sqrt (l  +  kron(m-l,2))*.  .  . 

sqrt ((n-1) “2-(abs (m-1) -1) ~2) *sqrt (n- 

-1+abs (m-1) ) 

*sqrt (n-l+abs (m-1) -2) ; 

59  g(n,m)=(m-l)/2*sqrt((2*(n-l)+l)/(2*(n-l)-l))*sqrt(n-l+ 

abs(m-l))*. . . 

sqrt (n-1- abs (m-1) ) ; 
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h(ii,m)  =  (m-l)  /4  /  abs  (m-1)  *sqrt  (  (2*(n-l)  +1)  /  (2*(n-l)  -1)  )  * 
sqrt((n-l) ~2-(abs(m-l)+l) "2)*. . . 

sqrt (n-l-abs (m-1) ) *sqrt (n-l-abs (m-1) -2) ; 


end 
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74 


if  (abs (m-1) ==1) 

d  (n , m) =0 ; 

g(n,m) =(m-l) /4/ abs (m-1) *sqrt ( (2*(n-l) +1) / (2*(n-l)  -1) ) * 
sqrt (n) *sqrt (n-2) *(n  +  l)  ; 

h(n,m) =(m-l) /4/ abs (m-1) *sqrt ( (2*(n-l) +1) / (2*(n-l)  -1) ) * 
sqrt (n-4) *  sqrt (n-3) *  .  .  . 
sqrt (n-2) *sqrt (n+1) ; 


end 

79  if  (abs (m-1) ==0) 

beta (n , m) =0 ; 

psi (n,m)=-(n  +  l) *sqrt ((n-1) *n/2)  ; 

84  elseif  (1  <  abs(m-l)  <  (n-1)) 

beta(n,m)  =  (n  +  l) /2*sqrt (1  +  A(abs (m-1)  ,  1) ) *sqrt (n-l  +  abs (m 
-l))*sqrt(n-l-abs(m-l)+l) ; 

psi (n,m)=-(n+l) /2*sqrt (n-l-abs (m-1) ) *sqrt (n-l+abs (m-1) 

+  1)  : 


89  end 

if  (abs(m-l)>  0) 
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mu (n,m)=-(m-l)/abs(m-l)*(n+l)/2*sqrt((2*(n-l)+l)/(2*(n 
-l)-l))*sqrt(l  +  A(abs(m-l)  ,1))*.  .  . 

94  sqrt(n-l  +  abs(m-l))*sqrt(ii-l  +  abs(m-l)-l); 
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eta(n,m)=-(m-l)/abs(m-l)*(n+l)/2*sqrt((2*(n-l)+l)/(2*( 
sqrt (n-l-abs (m-1) ) *sqrt (n-l-abs (m-1) -1) ; 


end 


end 


end 


Listing  A. 4: 

7o  This  f  unct  ion  comput  es  the  normal  field  potential  coefficients 
function  [J] =Normal (GM , AX , Nmax) 

GMS=0 . 3986005el5 ; 

4  AXS=6378137 . 0 ; 

JJ=0. 108262982131e-2; 

F1NV=298 . 257 ; 

FLTN=1 . 0/FlNV ; 

E2=FLTN* (2 . 0-FLTN) ; 

9 

J(1)=GMS/GM; 

J(3)  =-0 .484169650276e-3; 

J(5)  =  0.790314704521e-6; 

J(7)  =-0 . 168729437964e-8; 

14  J(9)  =  0 . 346071647263e-ll ; 

J(ll)=-0.265086254269e-14; 
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Listing  A. 5: 


I  7. 

7,  This  function  computes  the  Legendre  function  and  its7o 
7,  first  order  derivatives  using  recursive  formulae  7o 

7.  7. 
7.  INPUT  7. 
7„  phi  =  latitude  of  the  desired  point  °L 
7„  Nmax  =  maximum  desired  degree  °L 
7,  Mmax  =  maximum  desired  order  7o 

7.  7. 
7.  OUTPUT  7. 
7.  7. 

7«  pnm  =  Normalized  associated  Legendre  function  7o 
7„  dP  =  first  order  derivative  of  the  Normalized  7o 
7,  associated  Legendre  f  unct  ion  7o 

7.  7. 
7.  REFERENCES  7. 
7„  Hwang,  C.  and  M.J.  Lin,  (1998)  ,  Fast  Integration7o 
7,  of  low  orbiter  '  s  trajctory  perturbed  by  the  7o 
7,  Earth  non-sphericity.  Journal  of  Geodesy,  °L 
7.  vol  72:578-585  7. 
1  7. 

7„  Borre  ,  Kai  (2004)  ,  Geoid  Undulations  computed  by7o 
7«  EGM96  ,  report,  Aalborg  University  7o 

7.  7. 
I  7. 

7,  By  Mehdi  Eshagh  and  Ramin  Kiamehr  2006  °L 
7„  DivisionofGeodesy  7o 
7,  Royal  Institute  of  Technology  7o 
7,  Stckholm,  Sweden  7o 
7.  Email  :  eshaghOkth  .  se  7o 
7.  7. 
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function  [pnm  ,  dP]  =Pnm  (phi  ,  Nmax  ,  Mmax) 

nrow=Nmax+l ; npl=Mmax+l ; 
phii=phi*pi/180; 


x=sin (phii ) ; 
y=cos (phii ) ; 

pnm  (  1  ,  1)  =1 . 0  ; 

pnm (2,l)=sqrt (3.0)*x; 

pnm (2 ,2)=sqrt (3.0) *y; 

pnm (3,2)=sqrt (5.0)*pnm(2,2)*x; 

for  i =3 : np 1 
n=i - 1 ; 

pnm (i , i) =sqrt ( (n+0 . 5) /n) *pnm (i-l,i-l)*y; 

end 

k=npl -1 ; 
for  i =3 ; k 
n=i - 1 ; 

pnm (i+1 , i)=sqrt (2. 0*n+3) *pnm (i , i) *x ; 

end 

nml =npl -2 ; 
for  j  =1 : nml 
m=j -1 ; 
k=j+2; 
for  i=k : npl 
n=i -1 ; 

c=(2.0*n+1.0)/(n-m)/(n+m) ; 
cl  =  c*(2.0*n-1.0)  ; 
c2  =  c*(n+m-l) / (2*n-3) *(n-m-l)  ; 
c 1  =  sqrt (cl); 
c2  =  sqrt ( c2 )  ; 

pnm (i , j ) =cl*x*pnm (i-1 , j ) -c2*pnm (i-2 , j )  ; 

end 
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70 


end 


for  n=l ; Nmax -1 
for  m=l : n 

dP(n,m)=sqrt (((n-1) -(m-1) ) *((n-l) +(m-l) +1) *(1  +  A(m-1 ,m-l) ) *. .  . 
pnm (n,m+l) -(m-1) *. . . 

75  tan (phi ) *pnm (n  ,  m) )  ; 

end 

end 


Listing  A. 6: 

2  y,  Gradients  due  to  KC-10  Tanker  7, 

7,  (used  to  generate  white  noise)7o 

close  all 
7  clear  all 
clc 

7o  z  positive  "downward" 

12  G=6 . 67E-11 ; 

p  =  132;  7odensity  contrast  for  tanker 
Eotvos  =  lE-9;  7,  Eotvos  conversion 


a= [-3 . 05+25  3 . 05+25] ; 

7.  Point 

1, 

Permuted 

for 

7 

remaining 

points 

17  b=[-5-24.4  -5+24.4] ; 

7.  Point 

1, 

Permuted 

for 

7 

remaining 

points 

c=[-5-3.05  -5+3.05] ; 

7.  Point 

1, 

Permuted 

for 

7 

remaining 

points 

alt  =0 ; 
Ng  =  0; 
Eg  =  0; 

22 
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sumTxx  =0  ; 
f  or  Ctrl =1:2 

for  ctr 2  =  1  :  2 

for  ctr3  =  1 : 2 

27  sumTxx  =  sumTxx  +  ((-l)~(ctrl)*(-l)~(ctr2)*(-l)~(ctr3)) 

atarL2(((Ng-b(ctr2))  .*(alt-c(ctr3)))  ,((Eg-a(ctrl)).. 
.*((Eg-a(ctrl)) .~2+(Ng-b(ctr2)) .~2+(alt-c(ctr3))... 
.-2)  .-(.5)))  ; 

end 

end 

end 

32  Txx=G*p* sumTxx/Eotvos 


sumTxy  =0 ; 
f  or  Ctrl =1:2 
37  for  ctr 2  =  1 : 2 

for  ctr3  =  1 : 2 

sumTxy  =  sumTxy  +  ((-l)~(ctrl  +  ctr2  +  ctr3))  .*log(((alt-c(... 
ctr3 ) ) + ( ( Eg -a ( ctr 1 ) )  .  " 2+ ( Ng-b ( ctr 2 ) )  . ~ 2+ ( alt - c ( ctr3 
))  .-2)  .-(.5)))  ; 

end 


end 

42  end 


Txy=-G*p* sumTxy ; 


47  sumTxz=0; 

f  or  Ctrl =1:2 

for  ctr 2  =  1 : 2 

for  ctr3  =  1 : 2 
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sumTxz  =  sumTxz  +  ((-l)~(ctrl)*(-l)~(ctr2)*(-l)“(ctr3))  .* 
log(((Ng-b(ctr2))  +  ((Eg-a(ctrl))  . ~2+(Ng-b(ctr2) ) . .  . 
."2+(c(ctr3)-alt)  .~2)  .“(.5))); 

52  end 

end 

end 

Txz=-G*p* sumTxz/Eotvos ; 

57 


sumTyy=0 ; 
f  or  Ctrl =1:2 

for  ctr 2  =  1 : 2 

62  f  or  ctr3  =  1 : 2 

sumTyy  =  sumTyy  +  ((-l)~(ctrl)*(-l)“(ctr2)*(-l)“(ctr3))  .* 
atan2(((Eg-a(ctrl))  .*(alt-c(ctr3)))  ,((Ng-b(ctr2)). 
.*((Eg-a(ctrl)) .~2+(Ng-b(ctr2)) .~2+(c(ctr3)-alt).. 
.-2)  .-(.5)))  ; 

end 


end 


end 

67 

Tyy  =  G*p*  sumTyy /Eotvos 


sumTyz  =0 ; 

72  f or  Ctrl =1:2 

for  ctr 2  =  1 : 2 

for  ctr3  =  1 : 2 

sumTyz  =  sumTyz  +  ((-l)~(ctrl)*(-l)~(ctr2)*(-l)~(ctr3))  .  * 
log(((Eg-a(ctrl))  +  ((Eg-a(ctrl))  . ~2+(Ng-b(ctr2) ) . .  . 
."2+(c(ctr3)-alt) . ~2)  .“(.5))); 

end 

77  end 


end 
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Tyz=-G*p* sumTyz/Eotvos ; 


82  Tzz =- ( Tyy +Txx ) 


Listing  A. 7: 

7o  Terrain  Avoidance  Scenario  7. 

close  all 
clear  all 
clc 

8  vi=  [50  0  0]  ; 
speed=norm ( vi ) ; 

G=6 . 67E-11 ; 

p  =  2670;  7odensity  contrast  ground 
7o  p  =  1000;  7.  density  contrast  water  tower 

13  spacing=l ; 

N=  [0  :  spacing  :  1550]  ;  7o  setup  grid 
E=  [0  :  spacing  :  1550]  ;  "L  setup  grid 

Eotvos  =  lE-9;  7oEotvos  conversion 
18  update_rate =1 ; 

filter_cutoff  =  . 2 ; 

7o7o  25m  Cubic  Object  -  dimensions  permuted  for  each  obstacle  size, 
water  tower  gradients  calculated  similarly 
lengthl=25;  7«  object  base  (assumes  square) 

23  c=[0  25]  ;  7.  obstacle  height 

1 1 =0 : spacing / speed : ( round ( length (N) /2) -(lengthl/ spacing) /2-1) /  .  .  . 
speed;  7o  Truth  signal  time  array 
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t2=0:l/update_rate  :  (  round  (length(N)/2)-lengthl/2-l)/speed;  7,  time 
array  for  sensor  with  IHz  sampling  rate 

28  a=[1000  1000+ lengthl  ]  ;  7,  putting  obstacle  on  grid 

b=  [round  (  length  ( N) /2)  - lengthl /2  round  (  length  (N ) /2)  +  lengthl /2]  ;  7. 

putting  obstacle  on  grid 
7o  obstacle  height 

alt  =  15;  7o  10m  below  obstacle  top 
33  ipos=[0  length(N)/2  alt]; 

[Eg , Ng] =meshgr id (E , N) ; 

7.TXX 

38  sumTxx=0; 

for  Ctrl =1 : 2 
for  ctr 2  =  1 ; 2 

for  ctr3  =  1 : 2 

43  sumTxx  =  sumTxx  +  ((-l)~(ctrl)*(-l)“(ctr2)*(-l)“(ctr3))  .* 

atan2(((Ng-b(ctr2))  .*(alt-c(ctr3)))  ,((Eg-a(ctrl)). 
.*((Eg-a(ctrl)) ."2+(Ng-b(ctr2)) .~2+(alt-c(ctr3)).. 
.-2)  .-(.5)))  ; 

end 

end 

end 

48  Txx=G*p* sumTxx/Eotvos ; 

7.  Txy 
sumTxy  =0  ; 
f  or  Ctrl =1:2 
53  for  Ctr 2  =  1 : 2 

for  ctr3  =  1 : 2 
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sumTxy  =  sumTxy  +  ((-l)~(ctrl  +  ctr2  +  ctr3))  .*log(((alt-c(... 
ctr3))+((Eg-a(ctrl))  ."2+(Ng-b(ctr2))  .~2+(alt-c(ctr3 
))  .-2)  .-(.5)))  ; 

end 

end 

58  end 

Txy=-G*p* sumTxy ; 

Txy  (727 , 1051)  =-0 . 000001  ;  ”/«  Prevent  infinite  values 
Txy  (727 , 1001)  =0 . 000001  ;  Prevent  infinite  values 

63  Txy  (777 , 1051)  =0 . 000001  ;  Prevent  infinite  values 

Txy  (777 , 1001)  =-0 . 000001  ;  ”/«  Prevent  infinite  values 
Txy=Txy/Eotvos ; 

7.  Txz 

68  sumTxz=0; 

f  or  Ctrl =1:2 

for  ctr 2  =  1 : 2 

for  ctr3  =  1 : 2 

sumTxz  =  sumTxz  +  ((-l)~(ctrl)*(-l)~(ctr2)*(-l)“(ctr3)) 
log(((Ng-b(ctr2))  +  ((Eg-a(ctrl))  . ~2+(Ng-b(ctr2) )  .  .  . 
."2+(c(ctr3)-alt) .~2)  .“(.5))); 

73  end 

end 

end 

Txz=G*p* sumTxz / Eot VOS ; 

78 

“/.  Tyy 
sumTyy=0 ; 
f  or  Ctrl =1:2 

for  ctr 2  =  1 : 2 

83  f  or  ctr3  =  1 : 2 

sumTyy  =  sumTyy  +  ((-l)~(ctrl)*(-l)“(ctr2)*(-l)“(ctr3)) 

atan2(((Eg-a(ctrl))  .*(alt-c(ctr3)))  ,((Ng-b(ctr2)).. 
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.*((Eg-a(ctrl)) .~2+(Ng-b(ctr2)) .~2+(c(ctr3)-alt).. 
.-2)  .-(.5)))  ; 

end 

end 

end 

88 

Tyy=G*p* sumTyy /Eotvos ; 

7.  Tyz 
sumTyz  =0 ; 

93  f or  Ctrl =1:2 

for  ctr 2  =  1 : 2 

for  ctr3  =  1 : 2 

sumTyz  =  sumTyz  +  ((-l)~(ctrl)*(-l)~(ctr2)*(-l)~(ctr3)) 
log(((Eg-a(ctrl))  +  ((Eg-a(ctrl))  . ~2+(Ng-b(ctr2) ) . .  . 
."2+(c(ctr3)-alt) . ~2)  .“(.5))); 

end 

98  end 

end 

Tyz=G*p* sumTyz / Eot VOS ; 

103  7.  Tzz 

Tzz =- ( Tyy +Txx ) ; 

7,  Simulation 

sim  (  '  terrain_avoidance_sim  '  ,  0:. 01:30);  7orun  sim  for  30s 
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Figure  A. 2:  GGI  of  Justice  Simulink  Block  Diagram. 
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Figure  A. 3:  Terrain  Avoidance  Simulink  Block  Diagram. 
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Appendix  B.  Legendre  Funetion  Coeffieients  used  for  EGM96 

Gradient  Galeulations 


an^m  ^0?  ni  O5  1 


m,2 


an,m  =  - - ^ - '-\/rd  —  {m  —  l)^v^n  +  mGn  —  m  +  2,  2  <  m  <  n 


,  (n  +  m  +  l)(n  +  m  +  2) 

On,m  = - ^7 - — T - ,  m  =  0, 1 

2[m  +  1) 


+  3n  +  2 

On,m  = - ^7 - ,  2<m<n 


m,0 


Cn^rn  =  - - ^ —  (m  +  l)^^/n  —  m^/n  +  m  +  2,  m  =  0, 1 


c?i,m  =  ^  >/n^  —  (m  +  1)^  X  v^n  —  my^n  +  m  +  2,  2  <  m  <  n 


dn.m.  — 


7Ft  /  2fL  I  ]_ 

V  2n — +  ^m,2\/n^  -  {m-  +  mGn  +  m  -  2,  2  <  m  < 


Am  V  2n  —  1 


n 


TTT-  /  2/1  H"  1 

^n,m  =  + lVn-l(n  +  2),  m  =  l 


m  /2n  +  1  . -  . - 

5'n,m  ^  Y  Y  2^  -  1  ~  2<m<n 


Tfi  /  2/1  H"  1 

=  -  2V^  -  +  2,  m=l 


4m  V  2n  —  1 


m  /2n  +  1 

'  ^rj  m  — 


n,m  —  — Y  2 - 1  —  (m  +  1)^V^  ~  myn  —  m  —  2,  2  <  m  <  n 


dn,m  0)  ^  0 


n  +  2 


Pn,m  =  a/1  +  5m, -  m  +  1,  1  <  m  < 


n 


135 


,  ,  ,  n{n  + 1) 

ipn,m  =  -{n  +  2)\l - :: - ,  m  =  0 


Tl  2  . -  ! - 

Wn,ra  = - ~  m\/ H  +  771  +  1,  1<771<77 


^"n,m 


nm  +  2m  2n  +  1 


2m  V  277  —  1 


-  a/T+Xm 


77777  +  2m  2n  +  1  , -  , - - 

Vnm  = - \  - \/n  —  m^Jn  —  m  —  1,  0  <  777 

'  ’  2m  V  277  -  1  ^ 


5p,q  =  1,  p  =  q 


Sp,q  =  0,  p^q 


0  <  m  <  77 

<  77 
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Appendix  C.  Supplementary  Figures 
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Time  (s) 


(a)  Runs  1-6,  1000m,  GGI  1,  Rough  Terrain. 


Time  (s)  Time  (s) 

(b)  Runs  7-12,  1000m,  GGI  1,  Smooth  Terrain. 

Figure  C.l:  Runs  1-12,  1000m,  GGI  1. 
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(a)  Runs  13-18,  1000m,  GGI  2,  Rough  Terrain. 


Time  (s)  Time  (s) 

(b)  Runs  19-24,  1000m,  GGI  2,  Smooth  Terrain. 

Figure  C.2:  Runs  13-24,  1000m,  GGI  2. 
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(a)  Runs  25-30,  2500m,  GGI  1,  Rough  Terrain. 
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(b)  Runs  31-36,  2500m,  GGI  1,  Smooth  Terrain. 


Figure  C.3:  Runs  25-36,  2500m,  GGI  1. 
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(a)  Runs  37-42,  2500m,  GGI  2,  Rough  Terrain. 


Time  (s)  Time  (s) 


(b)  Runs  43-48,  2500m,  GGI  2,  Smooth  Terrain. 


Figure  C.4:  Runs  37-48,  2500m,  GGI  2. 
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(a)  Runs  49-54,  5000m,  GGI  1,  Rough  Terrain. 


Time  (s)  Time  (s) 

(b)  Runs  55-60,  5000m,  GGI  1,  Smooth  Terrain. 

Figure  C.5:  Runs  49-60,  5000m,  GGI  1. 
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(a)  Runs  61-66,  5000m,  GGI  2,  Rough  Terrain. 


Time  (s)  Time  (s) 

(b)  Runs  67-72,  5000m,  GGI  2,  Smooth  Terrain. 

Figure  C.6:  Runs  61-72,  5000m,  GGI  2. 
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(a)  Runs  73-78,  10000m,  GGI  1,  Rough  Terrain. 


r/-\A7\'7^-AAT^ 
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\f 

— 100  m/s  GGI 
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^■^0  20  40  60  80  100 


(b)  Runs  79-84,  10000m,  GGI  1,  Smooth  Terrain. 

Figure  C.7:  Runs  73-84,  10000m,  GGI  1. 
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(a)  Runs  85-90,  10000m,  GGI  2,  Rough  Terrain. 


Time  (s) 


(b)  Runs  91-96,  10000m,  GGI  2,  Smooth  Terrain. 

Figure  C.8:  Runs  85-96,  10000m,  GGI  2. 
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(a)  Runs  97-102,  20000m,  GGI  1,  Rough  Terrain. 
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(b)  Runs  103-108,  20000m,  GGI  1,  Smooth  Terrain. 


Figure  C.9:  Runs  97-108,  20000m,  GGI  1. 
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(a)  Runs  109-114,  20000m,  GGI  2,  Rough  Terrain. 
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(b)  Runs  115-120,  20000m,  GGI  2,  Smooth  Terrain. 

Figure  C.IO:  Runs  109-120,  20000m,  GGI  2. 
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(a)  Signal  time  rate  of  change,  noise  free,  rough  terrain. 
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(b)  Signal  time  rate  of  change,  noise  free,  smooth  terrain. 

Figure  C.ll:  Signal  Time  Rate  of  Change,  1000m. 
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(a)  Signal  time  rate  of  change,  noise  free,  rough  terrain. 


(b)  Signal  time  rate  of  change,  noise  free,  smooth  terrain. 

Figure  C.12:  Signal  Time  Rate  of  Change,  2500m. 
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(a)  Signal  time  rate  of  change,  noise  free,  rough  terrain. 


(b)  Signal  time  rate  of  change,  noise  free,  smooth  terrain. 

Figure  C.13:  Signal  Time  Rate  of  Change,  5000m. 
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signal  T^me  Rate  of_Change  (Eo_/s)  _  Signal  Time  Rate  of  Change  (Eo/s) 


(a)  Signal  time  rate  of  change,  noise  free,  rough  terrain. 


(b)  Signal  time  rate  of  change,  noise  free,  smooth  terrain. 


Figure  C.14:  Signal  Time  Rate  of  Change,  10000m. 
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(a)  Signal  time  rate  of  change,  noise  free,  rough  terrain. 


(b)  Signal  time  rate  of  change,  noise  free,  smooth  terrain. 

Figure  C.15:  Signal  Time  Rate  of  Change,  20000m. 
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(a)  GGI2  Tanker  Effect,  Rough  Terrain. 


(b)  GGI2  Tanker  Effect,  Smooth  Terrain. 

Figure  C.16:  Runs  123-124,  GGI2  Tanker  Effect,  5000m,  150m/s. 


Figure  C.17:  Runs  125-126,  GGIl  Tanker  Effect,  10000m,  150m/s. 
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(a)  GGI2  Tanker  Effect,  Rough  Terrain. 


Figure  C.18:  Runs  127-128,  GGI2  Tanker  Effect,  10000m,  150m/s. 
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Time  to  Terrain  Impact  (s) 

(a)  Signal. 


Time  to  Terrain  Impact  (s) 

(b)  Signal  rate  of  change. 

Figure  C.19:  Terrain  Avoidance  Scenario,  50  m/s,  25x25x25m  Block. 
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Time  to  Terrain  Impact  (s) 


(a)  Signal. 


Time  to  Terrain  Impact  (s) 

(b)  Signal  rate  of  change. 

Figure  C.20:  Terrain  Avoidance  Scenario,  50  m/s,  50x50x50m  Block. 
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Time  to  Terrain  Impact  (s) 


(a)  Signal. 


Time  to  Terrain  Impact  (s) 

(b)  Signal  rate  of  change. 

Figure  C.21:  Terrain  Avoidance  Scenario,  50  m/s,  lOOxlOOxlOOm  Block. 
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Time  to  Terrain  Impact  (s) 


(a)  Signal. 


Time  to  Terrain  Impact  (s) 

(b)  Signal  rate  of  change. 

Figure  C.22:  Terrain  Avoidance  Scenario,  50  m/s,  250x250x250m  Block. 
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Time  to  Terrain  Impact  (s) 


(a)  Signal. 


Time  to  Terrain  Impact  (s) 


(b)  Signal  rate  of  change. 

Figure  C.23:  Terrain  Avoidance  Scenario,  50  m/s,  500x500x500m  Block. 
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an  aircraft  and  combined  with  a  GGI,  form  the  basis  for  a  covert  navigation  system  using  a  map  matching  process. 
This  system  is  completely  passive  and  essentially  unjammable. 

To  determine  feasibility  of  this  method,  a  GGI  sensor  model  was  developed  to  investigate  signal  levels  at 
representative  flight  conditions.  Aircraft  trajectories  were  simulated  over  modeled  gravity  gradient  maps  to 
determine  the  utility  of  flying  modern  GGIs  in  the  roles  of  navigation  and  terrain  avoidance.  It  was  shown  that  the 
GGI  based  map-matching  navigation  system  can  likely  provide  a  marked  improvement  over  a  non-aided  INS  but  is 
limited  by  decreasing  gravity  gradient  strength  at  higher  altitudes,  particularly  over  smooth  terrain.  Additionally, 

GGI  output  rate  and  bandwidth  limitations,  along  with  the  inverse  nature  of  the  terrain  avoidance  problem, 
rendered  GGI  aided  terrain  avoidance  unfeasible  for  the  time  being. 
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